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Abstract 

Moving mesh methods (also called r-adaptive methods) are space-adaptive strategies used for 
the numerical simulation of time-dependent partial differential equations. These methods keep 
the total number of mesh points fixed during the simulation, but redistribute them over time to 
follow the areas where a higher mesh point density is required. There are a very limited number 
of moving mesh methods designed for solving field-theoretic partial differential equations, and 
the numerical analysis of the resulting schemes is challenging. In this paper we present two 
ways to construct r-adaptive variational and multisymplectic integrators for (l+l)-dimensional 
Lagrangian field theories. The first method uses a variational discretization of the physical 
equations and the mesh equations are then coupled in a way typical of the existing r-adaptive 
schemes. The second method treats the mesh points as pseudo-particles and incorporates their 
dynamics directly into the variational principle. A user-specified adaptation strategy is then 
enforced through Lagrange multipliers as a constraint on the dynamics of both the physical field 
and the mesh points. We discuss the advantages and limitations of our methods. Numerical 
results for the Sine-Gordon equation are also presented. 

1 Introduction 

The purpose of this work is to design, analyze and implement variational and multisymplectic 
integrators for Lagrangian partial differential equations with space-adaptive meshes. In this paper 
we combine geometric numerical integration and r-adaptive methods for the numerical solution of 
PDEs. We show that these two fields are compatible, mostly due to the fact that in r-adaptation 
the number of mesh points remains constant and we can treat them as additional pseudo-particles 
whose dynamics is coupled to the dynamics of the physical field of interest. 

Geometric (or structure-preserving) integrators are numerical methods that preserve geometric 
properties of the flow of a differential equation (see |20|). This encompasses symplectic integrators 
for Hamiltonian systems, variational integrators for Lagrangian systems, and numerical methods on 
manifolds, including Lie group methods and integrators for constrained mechanical systems. Geo- 
metric integrators proved to be extremely useful for numerical computations in astronomy, molec- 
ular dynamics, mechanics and theoretical physics. The main motivation for developing structure- 
preserving algorithms lies in the fact that they show excellent numerical behavior, especially for 
long-time integration of equations possessing geometric properties. 
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An important class of structure-preserving integrators are variational integrators for Lagrangian 
systems ([20]. |34|). This type of integrators is based on discrete variational principles. The vari- 
ational approach provides a unified framework for the analysis of many symplectic algorithms and 
is characterized by a natural treatment of the discrete Noether theorem, as well as forced, dis- 
sipative and constrained systems. Variational integrators were first introduced in the context of 
finite-dimensional mechanical systems, but later Marsden & Patrick & Shkoller |31j generalized this 
idea to field theories. Variational integrators have since then been successfully applied in many 
computations, for example in elasticity (|30|). electrodynamics (|48j) or fluid dynamics ([38]). The 
existing variational integrators so far have been developed on static, mostly uniform spatial meshes. 
The main goal of this paper is to design and analyze variational integrators that allow for the use 
of space-adaptive meshes. 

Adaptive meshes used for the numerical solution of partial differential equations fall into three 
main categories: /i-adaptive, p-adaptive and r-adaptive. i?-adaptive methods, which are also known 
as moving mesh methods ([7], |25|). keep the total number of mesh points fixed during the simulation, 
but relocate them over time. These methods are designed to minimize the error of the computations 
by optimally distributing the mesh points, contrasting with /i-adaptive methods for which the 
accuracy of the computations is obtained via insertion and deletion of mesh points. Moving mesh 
methods are a large and interesting research field of applied mathematics, and their role in modern 
computational modeling is growing. Despite the increasing interest in these methods in recent 
years, they are still in a relatively early stage of their development compared to the more matured 
/i-adaptive methods. 

Overview 

There are three logical steps to r-adaptation: 

• Discretization of the physical PDE 

• Mesh adaptation strategy 

• Coupling the mesh equations to the physical equations 

The key ideas of this paper regard the first and the last step. Following the general spirit of vari- 
ational integrators, we discretize the underlying action functional rather than the PDE itself, and 
then derive the discrete equations of motion. We base our adaptation strategies on the equidis- 
tribution principle and the resulting moving mesh partial differential equations (MMPDEs). We 
interpret MMPDEs as constraints, which allows us to consider novel ways of coupling them to the 
physical equations. Note that we will restrict our explanations to one time and one space dimension 
for the sake of simplicity. 

Let us consider a (l+l)-dimensional scalar field theory with the action functional 



where <j) : [0,V maai ] x [0,T max ] — > R is the field and C : M x M. x R — > M its Lagrangian density. 
For simplicity, we assume the following fixed boundary conditions 
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<Ko,t) = 0l 

4>(X ma x,t) = 4>R 



(1.2) 



In order to further consider moving meshes let us perform a change of variables X = X(x, t) such 
that for all t the map X(.,t) : [0, X max ] — > [0,X max ] is a 'diffeomorphism' — more precisely, we 
only require that X(.,t) is a homeomorphism such that both X(.,t) and X(.,t)~ 1 are piecewise C 1 . 
In the context of mesh adaptation the map X(x, t) is going to represent the spatial position at time 
t of the mesh point labeled by x. Define <p(x,t) = cp(X(x, t), t). Then the partial derivatives of (j) 
are (j>x(X(x,t),t) = ip x /X x and cj>t(X(x,t),t) = <ft — tp x Xt/X x . Plugging these equations in 
we get 

rTm ax rX max fin in X* \ ~ 

S[<f>] = / cU^,<p t -^-t)X x dxdt=:S[cp},S[<p,X} (1.3) 

JO JO V A z ^x 7 

where the last equality defines two modified, or 'reparametrized', action functionals. For the first 
one, S is considered as a functional of ip only, whereas in the second one we also treat it as a 
functional of X. This leads to two different approaches to mesh adaptation, which we dub the 
control-theoretic strategy and the Lagrange multiplier strategy, respectively. The 'reparametrized' 
field theories defined by S[ip] and S[ip, X] are both intrinsically covariant; however, it is convenient 
for computational purposes to work with a space-time split and formulate the field dynamics as an 
initial value problem. 

Outline 

This paper is organized as follows. In Section[2]and Section|3]we take the view of infinite dimensional 
manifolds of fields as configuration spaces, and develop the control-theoretic and Lagrange multiplier 
strategies in that setting. It allows us to discretize our system in space first and consider time 
discretization later on. It is clear from our exposition that the resulting integrators are variational. 
In Section U we show how similar integrators can be constructed using the covariant formalism of 
multisymplectic field theory. We also show how the integrators from the previous sections can be 
interpreted as multisymplectic. In Section [5] we apply our integrators to the Sine-Gordon equation 
and we present our numerical results. We summarize our work in Section [6] and discuss several 
directions in which it can be extended. 

2 Control-theoretic approach to r-adaptation 

At first glance, it appears that the simplest and most straightforward way to construct an r-adaptive 
variational integrator would be to discretize the physical system in a similar manner to the general 
approach to variational integration, i.e. discretize the underlying variational principle, and then 
derive the mesh equations and couple them to the physical equations in a way typical of the existing 
r-adaptive algorithms. We explore this idea in this section and show that it indeed leads to space 
adaptive integrators that are variational in nature. However, we also show that those integrators 
do not exhibit the behavior expected of geometric integrators, such as good energy conservation. 
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2.1 Reparametrized Lagrangian 

For the moment let us assume that X(x,t) is a known function. We denote by £(X,t) the function 
such that f(.,t) = X(.,t)-\ that is £(X(x,t),t) = x0. We thus have = S[(p(£(X,t),t)]. 

Proposition 2.1. Extremizing S[<j>] with respect to (j) is equivalent to extremizing S[<p] with respect 
to if. 

Proof. The variational derivatives of S and S are related by the formula 

6S[<p] ■ 6<p(x, t) = 6S[<p{Z{X, t),t)} ■ 6ip(£{X, t),t). (2.1) 

Suppose (f)(X,t) extremizes S [</)], i.e. SS[<f>] -5(f) = for all variations Sep. Given the function X(x,t), 
define ip(x,t) = cp(X(x, t), t). Then by the formula above we have <W[y] = 0, so (f extremizes S. 
Conversely, suppose (f(x, t) extremizes S, that is 5S[ip]-5(p = for all variations 5ip. Since we assume 
X(.,t) is a homeomorphism, we can define cj)(X,t) = <p(£(X,t),t). Note that an arbitrary variation 
5(p(X,t) induces the variation 5(p(x,t) = 5(fr(X(x,t),t). Then we have SS[4>] ■ 5(f) = SS[(p] ■ 5<p = 
for all variations 5(j>, so <p(X,t) extremizes S[<fi]. 

□ 

The corresponding instantaneous Lagrangian L : Q x W x M — ► R is 

L[tp,(p t ,t] = / £((p,tp x ,tp t ,t)dx (2.2) 
Jo 

with the Lagrangian density 

C((p,(p x ,ip t ,x,t) = c(p, Y-,ft - ^y-^X x . (2.3) 

The function spaces Q and W must be chosen appropriately for the problem at hand, so that (|2.2p 
makes sense. For instance, for a free field we will have Q = H 1 ([0, X rnax ]) and W = L 2 ([0, X max ]). 
Since X(x,t) is a function of t, we are looking at a time-dependent system. Even though the energy 
associated with (|2.2p is not conserved, the energy of the original theory associated with (II. ip 




is conserved. To see this, note that if 4>{X,t) extremizes S[(j)} then dE/dt = (computed from 
(12. 4p ). Trivially, this means that dE/dt = when formula (12. 5|) is invoked as well. Moreover, as we 
have noted earlier, cj)(X,t) extremizes S[cp] iff <p(x,t) extremizes S[(f\. This means that the energy 
(|2,5p is constant on solutions of the reparametrized theory. 

1 We allow a little abuse of notation here: X denotes both the argument of £ and the change of variables X(x,t). 
If we wanted to be more precise, we would write X = h(x,t). 
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2.2 Spatial Finite Element discretization 



We begin with a discretization of the spatial dimension only, thus turning the original infinite- 
dimensional problem into a time-continuous finite-dimensional Lagrangian system. Let Ax = 
X ma x/{N + 1) and define the reference uniform mesh Xj = i ■ Ax for i = 0, 1, ...,N + 1, and 
the corresponding piecewise linear finite elements 



m{x) = { -^i, ifxi<x<x m , (2.6) 



0, otherwise. 
We now restrict X(x, t) to be of the form 



N+l 

X{x,t)=Y J X i {t)i li {x) (2.7) 

i=0 

with Xo(t) = 0) Xjy + i(t) = X max and arbitrary Xi(t), i = 1,2, ...,N as long as X(.,t) is a 
homeomorphism for all t. In our context of numerical computations, the functions Xi(t) represent 
the current position of the i th mesh point. Define the finite element spaces 



Qn = W N = span (770, ...,r] N +i) ( 2 - 8 ) 

and assume that Qn C Q, Wn C W. Let us denote a generic element of Qn by ip and a generic 
element of Wn by 99. We have the decompositions 

N+l N+l 

^(x) = y^ x ) = Y y^ x )- ( 2 - 9 ) 

i=0 i=0 

The numbers (j/j, y,) thus form natural (global) coordinates on Qn x Wn- We can now approximate 
the dynamics of system (I2.2p in the finite-dimensional space Qn x Wat- Let us consider the restriction 
Ln = -^IqjvxVKjvxR °f the Lagrangian (|2.2p to Qn x Wn X M. In the chosen coordinates we have 

N+l N+l 

L N {yo,.-.,yN+i,m,---,yN+i,t) = L^^2y i r] i (x),J2y i ri i (x),t . (2.10) 

i=0 i=0 

Note that, given the boundary conditions (jl.2p . yo, yN+i, ilo, an d i/N+i are fixed. We will thus no 
longer write them as arguments of Ln- 

The advantage of using a finite element discretization lies in the fact that the symplectic structure 
induced on Qn x Wn by Ln is strictly a restriction (i.e., a pull-back) of the (pre-)symplectic 
structured on Q x W. This establishes a direct link between symplectic integration of the finite- 
dimensional mechanical system (Qn x Wn, Ln) and the infinite-dimensional field theory (Q x W, L) 



2 In most cases the symplectic structure of (Q x W, L) is only weakly-nondegenerate; see [16] 
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2.3 DAE formulation and time integration 

We are now going to consider time integration of the Lagrangian system (Qn X Wn,Ln). If the 
functions Xi(t) are known, then one can perform variational integration in the standard way, that 
is, define the discrete Lagrangian L,j:Rx Qn X R x Qn — >• R and solve the corresponding discrete 
Euler-Lagrange equations (see |34j . |20|). Let i n = n • At for n = 0,1,2,... be an increasing 
sequence of times and {y , y , . . .} the corresponding discrete path of the system in Q^r. The 
discrete Lagrangian L d is an approximation to the exact discrete Lagrangian L^, such that 

/•tn+1 

L d (t n ,y n ,t n+l ,y n+l ) « L^(t n ,y n ,t n+l ,y n+l ) = / L N (y(t),y(t),t) dt, (2.11) 

•/tn 

where y n = (y™, y^v), y n+1 = (y™ +1 , y^ +1 ) and y(i) is the solution of the Euler-Lagrange 
equations corresponding to Ljy with the boundary values y(t n ) = y n , y(t n+ i) = y n+1 . Depending on 
the quadrature we use to approximate the integral in (I2.1ip , we obtain different types of variational 
integrators. As will be discussed below, in r-adaptation one has to deal with stiff differential 
equations or differential-algebraic equations, therefore higher order integration in time is required. 
We are going to employ variational partitioned Runge-Kutta methods. An s-stage Runge Kutta 
method is constructed by choosing 

s 

L d (t n , y n ,t n+1 ,y n+1 ) = (t n+1 - t n ) biL N (Yi, Y h U), (2.12) 

i=i 

where t{ = t n + Cj(t n +i — t n ), the right-hand side is extremized under the constraint y n+1 = 
y n + (tn+i — tn) Yli=i biYi, and the internal stage variables Yi, Y{ are related by Yi = y n + (t n +i — 
tn) Ylj=i a ijYj- R can be shown that the variational integrator with the discrete Lagrangian (|2.12p 
is equivalent to an appropriately chosen symplectic partitioned Runge-Kutta method applied to 
the Hamiltonian system corresponding to Ljy (see [34] . |2U|). With this in mind we turn our semi- 
discrete Lagrangian system (Qn x Wn,Ln) into the Hamiltonian system (Qn x W^,H^) via the 
standard Legendre transform 

N 

Hn(vi, ...,yN,Pi, -,Pn;X 1 , ...,X n ,Xi, ...,X N ) = ^Piiji - L N (yi,...,y N ,y 1 , ...,y N ,t), (2.13) 

i=l 

where pi = dL^/dyi and we explicitly stated the dependence on the positions Xj and velocities Xj 
of the mesh points. The Hamiltonian equations take the forrrH 

m = ^(y, P ;X(t),X(t)), (2.14) 
Pi = -^§f(y,p;X(t),x(t)). 

3 It is computationally more convenient to directly integrate the implicit Hamiltonian system pt = OLn /dyt, 
pi = OLn /dyi, but as long as system (II. 1|) is at least weakly-nondegenerate there is no theoretical issue with passing 
to the Hamiltonian formulation, which we do for the clarity of our exposition. 
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Suppose that the functions Xi(t) are C 1 and Hjy is smooth as a function of the y^s, PiS, X^s and 
Xj's (note that these assumptions are used for simplicity, and can be easily relaxed if necessary, 
depending on the regularity of the considered Lagrangian system). Then the assumptions of Picard's 
theorem are satisfied and there exists a unique C 1 flow Ft 0: t = (F]? t , Ff t ) : Qn x Wjy — > Qn x W^- 
for (|2.14|) . This flow is symplectic. 

However, in practice we do not know the Xj's and we in fact would like to be able to adjust them 
'on the fly', based on the current behavior of the system. We are going to do that by introducing 
additional constraint functions gi(yi, yjy, X%, Xpf) and demanding that the conditions gi = 
be satisfied at all timea^l The choice of these functions will be discussed in Section 12.41 This leads 
to the following differential-algebraic system of index 1 (see [BJ, |22j . |23j) 



m = ^(y,p-,X,x), (2.15) 

Pi = — Q^~{ y > p ' x > x )> 

o = g t (y,x), 
yi(to) = yf\ 

Pi{h) =Pi 

for i = 1,...,N. Note that an initial condition for X is fixed by the constraints. This system is of 
index 1, because one has to differentiate the algebraic equations with respect to time once in order 
to reduce it to an implicit ODE system. In fact, the implicit system will take the form 



9H N , 

v = — n — [y,p;X,x 



y = ^-[y.i>:X.X). (2.1(5) 

91 

dy 

= ^(y,X)y+^(y,X)X, 

y(t ) = y {0 \ 
P(t )=p {0) , 
X(t ) = x(°\ 

where is a vector of arbitrary initial condition for the JQ's. Suppose again that Hn is a smooth 
function of y, p, X and X. Futhermore, suppose that g is a C 1 function of y, X, and ~ 8y ^ s 
invertible with its inverse bounded in a neighborhood of the exact solution^ Then, by the Implicit 
Function Theorem equations (I2.16P can be solved explicitly for y, p, X and the resulting explicit 
ODE system will satisfy the assumptions of Picard's theorem. Let (y(t),p(t),X(t)) be the unique 
C 1 solution to this ODE system (and hence to (|2.16p ). We have the trivial result 

4 In the context of Control Theory the constraints gi — are called strict static state feedback. See |37] . 
5 Again, these assumptions can be relaxed if necessary. 
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Proposition 2.2. If g(y^,X^) = 0, then (y(t),p(t), X(t)) is a solution to [2J5\) H 

In practice we would like to integrate system (|2.15p . A question arises in what sense is this system 
symplectic and in what sense a numerical integration scheme for this system can be regarded as 
variational. Let us address these issues. 

Proposition 2.3. Let (y(t),p(t),X(t)) be a solution to A2.15\) and use this X(t) to form the Hamil- 
tonian system \2.1Jt ). Then we have that 



t ,t\y if )■> yy) — * t , 

and 



9{Fl t (yWJV),X(t))=Q, 
where F to j(y,p) is the symplectic flow for \2.1J$. 



Proof. Note that the first two equations of (|2.15p are the same as (|2.14p . therefore (y(t),p(t)) 
trivially satisfies (|2.14|) with the initial conditions y(to) = y^ and p(to) = p(°\ Since the flow map 
F^t is unique, we must have y(t) = Ff Q>t {y^\p^) and p(t) = F? o t (y(°\p^). Then we also must 

have that g(^F^ t {y^\p^), X{t)^ = 0, that is, the constraints are satisfied along one particular 

integral curve of (|2.14p that passes through (2/ (0) ,p> (0) ) at t Q . 

□ 

Suppose we now would like to find a numerical approximation of the solution to (|2,14p using 
an s-stage partitioned Runge-Kutta method with coefficients Ojj, a^, hi, Ci (|21{. |2U|). The 
numerical scheme will take the form 



P { = — ^ ( Y\ P l ;X(t n + CiAt),X(t n + Cl At) 
dy 



Y* = y n + AtY,a lJ Yi, 

s 

P l =p n + AtY / a lJ P j , 

i=i 

s 

y n+1 = y n + AtY,b t Y\ 



p n+l = p 



^ + AtY,hP\ 



i=i 



6 Note that there might be other solutions, as for any given y(°' there might be more than one X' ' that solves 
the constraint equations. 



s 



where Y l , Y'\ P l , P l are the internal stages and At is the integration timestep. Let us apply the 
same partitioned Runge-Kutta method to (I2.15p , In order to compute the internal stages Q l , Q l of 
the X variable we use the state-space form approach, that is, we demand that the constraints and 
their time derivatives be satisfied (see |22|). The new step value X n+l is computed by solving the 
constraints as well. The resulting numerical scheme is thus 



yi = ^EHi Y \P l ;Q\Q t ), (2.18) 



dp 

P i = -^-(y^P^Q^Q' 1 
ay v 

s 

Y i = y n + AtJ2aijY 3 , 
i=i 

s 

P* =p n + AtY j a il P 3 , 

3=1 

= g(Y i ,Q i ), 

= S (yi ' gi) ^ + ll (yi ' Ql) ^ 



y n+1 = y n + AtJ2biY i , 

i=l 



p n + 1 =p n + AtJ2b i P i , 

i=l 
n+1 vn+l\ 



= g(y n+ \X 
We have the following trivial observation. 

Proposition 2.4. If X(t) is defined to be a C 1 interpolation of the internal stages Q l , Q l at times 
t n + CiAt (that is, if the values X(t n + CiAt), X{t n + CiAt) coincide with Q l , Q l ), then the schemes 
112.17]) and 112.18]) give the same numerical approximations y n , p n to the exact solution y(t), p(t). 

Intuitively Proposition 12.41 states that we can apply a symplectic partitioned Runge-Kutta 
method to the DAE system (|2.15p . which solves both for X{t) and (y(t),p(t)), and the result will be 
the same as if we performed a symplectic integration of the Hamiltonian system (|2,14p for (y(t),p(t)) 
with a known X(t). 



2.4 Moving mesh partial differential equations 

The concept of equidistribution is the most popular paradigm of r-adaptation (see [7|, |25|). Given 
a continuous mesh density function p(X), the equidistribution principle seeks to find a mesh = 
Xq < X\ < ... < Xjv+i = X max such that the following holds 

/ p{X)dX= p{X)dX = ...= l p(X)dX, (2.19) 

JO JX X JX N 
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that is, the quantity represented by the density function is equidistributed among all cells. In the 
continuous setting we will say that the reparametrization X = X(x) equidistributes p(X) if 



X(x) x 



p(X) dX = a, (2.20) 



o X max 



where a = J max p(X) dX is the total amount of the equidistributed quantity. Differentiate this 
equation with respect to x to obtain 

r)X 1 

p(X(x))— = —a. (2.21) 

OX -A-max 

It is still a global condition in the sense that a has to be known. For computational purposes it is 
convenient to differentiate this relation again and consider the following partial differential equation 

d / . . . , dX \ 

with the boundary conditions X(0) = 0, X(X max ) = X max . The choice of the mesh density function 
p(X) is typically problem-dependent and the subject of much research. A popular example is the 
generalized solution arclength given by 



^ 2 (l)V i+ ° 2 (t) 2 (2 - 23) 

It is often used to construct meshes that can follow moving fronts with locally high gradients ([7], 
|25|). With this choice, equation (|2,22p is equivalent to 

a 2 f x fxx + X X X XX = 0, (2.24) 

assuming X x > 0, which we demand anyway. A finite difference discretization on the mesh X{ = i-Ax 
gives us the set of contraints 



9i(yi, -,yN,Xi, -,x N ) = 

a 2 (y i+1 - yi ) 2 + (X i+ i - Xi) 2 - a 2 ( yi - y^) 2 - (X, - X^) 2 = 0, (2.25) 
with the previously defined yi's and AYs. This set of constraints can be used in (|2.15p . 

2.5 Example 

To illustrate these ideas let us consider the Lagrangian density 

C{4>,4> x ,4>t) = \<i> 2 -W{<j>x). (2.26) 
The reparametrized Lagrangian (|2.2I) takes the form 

i \ 2 / ,n \ 1 

dx. (2.27) 



L[(p,ip t ,t] 



o 



\x x (<& - ^x t ) 2 - w( ^ ) x x 
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Let N = 1 and 4>l = 4>R = 0- Then 



tp(x,t) = y 1 (t)rj 1 (x), X(x,t)=X 1 (t)rj 1 (x) + X maxm (x). (2.28) 

The semi-discrete Lagrangian is 

f ( ■ A X ^)(- ^1 v i X max -X x (t) ( . y± ■ ^ 2 

Ln{VU * } =— [ yi ~ xW) Xl{t) ) + 6 [ yi + X max -X l{ t) X ^ 



The Legendre transform gives p\ = dL^ /diji = X max yi/3, hence the semi-discrete Hamiltonian is 



£V i v v \ 3 2 1 X max X\ 2 
H N (y 1 ,pi;X 1 ,X 1 ) =— Pl - -^ryr; tttVi 



The corresponding DAE system is 



m = it—pi, (2-31) 



X m ao- 



1 X max Xf ur>fyi\ , txt'/ W 



3 Xi (X max — X\ ) V Xi / V X max — X\ 

= 5l (y 1 ,X 1 ). 

This system is to be solved for the unknown functions yi(t), pi(t) and X\(t). It is of index 1, because 
we have three unknown functions and only two differential equations — the algebraic equation has 
to be differentiated once in order to obtain a missing ODE. 



2.6 Backward error analysis 

The true power of symplectic integration of Hamiltonian equations is revealed through backward 
error analysis: it can be shown that a symplectic integrator for a Hamiltonian system with the 
Hamiltonian H(q,p) defines the exact flow for a nearby Hamiltonian system, whose Hamiltonian 
can be expressed as the asymptotic series 



M> (q,p) = H{q,p) + AtH 2 (q,p) + At 2 H 3 (q,p) + ... (2.32) 

Owing to this fact, under some additional assumptions symplectic numerical schemes nearly conserve 
the original Hamiltonian H(q,p) over exponentially long time intervals. See [20] for details. 

Let us briefly review the results of backward error analysis for the integrator (12, 18ft . Suppose 
g(y,X) satisfies the assumptions of the Implicit Function Theorem. Then, at least locally, we can 
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solve the constraint X = h(y). The Hamiltonian DAE system f|2. 15j) can be then written as the 
following (implicit) ODE system for y and p 



Since we used the state-space formulation, the numerical scheme (|2. 18|) is equivalent to applying the 
same partitioned Runge-Kutta method to (|2,33p . that is, we have Q l = h(Y l ) and Q l = h'(Y t )Y' 1 . 
We computed the corresponding modified equation for several symplectic methods, namely Gauss 
and Lobatto IIIA-IIIB quadratures. Unfortunately, none of the quadratures resulted in a form akin 
to ([2.33p for some modified Hamiltonian function related to Hn by a series similar to (|2.32|) . 
This hints at the fact that we should not expect this integrator to show excellent energy conservation 
over long integration times. One could also consider the implicit ODE system (12,16|) . which has an 
obvious triple partitioned structure, and apply a different Runge-Kutta method to each variable y, 
p and X. Although we did not pursue this idea further, it seems unlikely it would bring a desirable 
result. 

We therefore conclude that the control-theoretic strategy, while yielding a perfectly legitimate 
numerical method, does not take the full advantage of the underlying geometric structures. Let us 
point out that, while we used a variational discretization of the governing physical PDE, the mesh 
equations were coupled in a manner that is typical of the existing r-adaptive methods (see [7], [25J). 
We now turn our attention to a second approach, which offers a novel way of coupling the mesh 
equations to the physical equations. 

3 Lagrange multiplier approach to r-adaptation 

As we saw in Section [2l discretization of the variational principle alone is not sufficient if we would 
like to accurately capture the geometric properties of the physical system described by (11. ip . In 
this section we propose a new technique of coupling the mesh equations to the physical equations. 
Our idea is based on the observation that in r-adaptation the number of mesh points is constant, 
therefore we can treat them as pseudo-particles, and we can incorporate their dynamics into the 
variational principle. We show that this strategy results in integrators that much better preserve 
the energy of the considered system. 

3.1 Reparametrized Lagrangian 

In this approach, we treat X(x, t) as an independent field, that is, another degree of freedom, and 
we are going to treat the 'modified' action (|1.3p as a functional of both 93 and X: S = S[cp,X], 
For the purpose of the derivations below, we assume that (p(.,t) and X(.,t) are continuous and 
piecewise C . One could consider the closure of this space in the topology of either Hilbert or 
Banach space of sufficiently integrable functions and interpret differentiation in a sufficiently weak 
sense, but this functional-analytic aspect is of little importance for the developments in this section. 
We refer the interested reader to |12j and |13| . As in Section [2.1[ let £(X, t) be the function such 




(2.33) 
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that = X(.,t)~ l , that is £(X(x,t),t) = x. Then S[<p,X] = S[<p(£(X,t),t)]. We begin with 

two propositions and one corollary which will be important for the rest of our exposition. 

Proposition 3.1. Extremizing S[<j)] with respect to <p is equivalent to extremizing S[cp,X] with 
respect to both tp and X. 

Proof. The variational derivatives of S and S are related by the formula 

SxS[(p, X] ■ 8ip(x, t) = 6S[<p(t(X, t),t)] ■ 8<p(t(X, t),t), (3.1) 

5 2 S[v, X] • SX(x, t) = 6S[<p(t(X, *),*)]■(- xlft&Xt) 5Xi * {X ' t} ' ' 

where 81 and 82 denote differentiation with respect to the first and second argument, respectively. 
Suppose 4>(X, t) extremizes S[<j)], i.e. 8S[<p] -Sep = for all variations 8(f). Choose an arbitrary X(x, t), 
such that X(.,t) is a (sufficiently smooth) homeomorphism and define <p(x,t) = <j)(X(x, t),t). Then 
by the formula above we have 8\S[<p,X] = and 82S[<p,X] = 0, so the pair {ip,X) extremizes 
S. Conversely, suppose the pair ((p,X) extremizes S, that is 8\S[ip,X] ■ Sip = and 52S[<p,X] ■ 
5X = for all variations 5<p and 5X. Since we assume X(.,t) is a homeomorphism, we can 
define (j)(X,t) = tp(£(X,t),t). Note that an arbitrary variation S(j>(X,t) induces the variation 
Scp(x,t) = 5<j)(X(x,t),t). Then we have 8S[<f>] ■ 8(f) = 8iS[<p, X] ■ 8<p = for all variations 8(f), so 
(j>(X,t) extremizes S[(f>]. 

□ 

Proposition 3.2. The equation 8 2 S[(f,X] = is implied by the equation 8iS[<p,X] = 0. 

Proof. As we saw in the proof of Proposition 13.11 the condition 8iS[ip,X] ■ Sep = implies 8S = 0. 
By (|3.ip . this in turn implies 82S[(p,X] ■ SX = for all 8X. Note that this argument cannot be 
reversed: 82§[tp, X] • 8X = does not imply 8S = when (p x = 0. 

□ 

Corollary 3.3. The field theory described by S[<p,X] is degenerate and the solutions to the Euler- 
Lagrange equations are not unique. 

3.2 Spatial Finite Element discretization 

The Lagrangian of the 'reparametrized' theory L: QxGxWxZ — > R 

L[tp,X,tp t ,X t ] = C(ip,-^-,ip t ^^)X x dx (3.2) 

has the same form as (I2.2p (we only treat it as a functional of X and X% as well), where Q, G, 
W and Z are spaces of continuous and piecewise C 1 functions, as mentioned before. We again let 
Ax = X max /(N + 1) and define the uniform mesh X{ = i ■ Ax for i = 0, 1, ...,N + 1. Define the 
finite element spaces 

Qn = Gn = W N = Z N = span(r/ ,...,r?jv+i), (3.3) 
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where we used the finite elements (|2.6p . We have Qn C Q, Gm C G, C W, Zjv C Z. In 
addition to (12.91) we also consider 



N+l 

X(x) = Xirhi? 

i=0 



7V+1 

X{x) = Y,XiVi{x)- 

i=Q 



(3.4) 



The numbers (yj, Xi,yi, Xj) thus form natural (global) coordinates on Q^r x Gat x VFjv x ^jv- We 
again consider the restricted Lagrangian = L\q nX g n xW n xZ n - In the chosen coordinates 



L N (y 1 ,...,y N ,X 1 , ...,X N ,yi, ...,y N ,Xi, ...,X N ) = L <p(x),X(x),<p(x),X(x) 



(3.5) 



where (p(x), X(x), (f{x), X(x) are defined by (|2.9p and f|3.4j) . Once again, we refrain from writing 
?/0) Vn+i, Vo, VN+1-, Xq, Xn+i, Xq and -Xjv+i as arguments of Ljy hi the remainder of this section, 
as those are not actual degrees of freedom. 

3.3 Invertibility of the Legendre Transform 

For simplicity, let us restrict our considerations to Lagrangian densities of the form 

1 



c(4>, <f> x , <t>t) = ^<Pt - R{<i>x, <t>)- 



(3.6) 



We chose a kinetic term that is most common in applications. The corresponding 'reparametrized' 
Lagrangian is 



L[<p,X, <p t ,X t ] 



Xmax 



-X x (^p t - j^X t I d.r 



(3.7) 



where we kept only the terms that involve the velocities (ft and Xt- The semi-discrete Lagrangian 
becomes 



N 



Xi+i — Xi 



i=o 



y-i+i - Vi 

Xi+i — Xi 



x L 



+ Vi+i 



Vi+i - Vi 
Xi+i — Xi 

_ Vi+i - Vi 
Xi+i — x t 



X i Vi+l 



Vi+i - Vi 
Xi+i — Xi 



X 



i+1 



Xi+i 



(3i 



Let us define the conjugate momenta via the Legendre Transform 



dL 



Pi 



N 



dm 



Si 



dL 



N 



dXi 



i = 1,2,..., N. 



(3.9) 



This can be written as 
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/ PI \ 

Si 



PN 

\S N J 



M N (y,X) 



( VI \ 

Xi 



VN 

\x N ) 



(3.10) 



where the 2N x 2N mass matrix Mj^(y,X) has the following block tridiagonal structure 

\ 



M N (y,X) 



( Ax Bi 
B\ A2 B2 

B 2 A 3 B% 



Bn-i 
Bn-i An J 



(3.11) 



with the 2x2 blocks 



where 

<5j = Xi+i - Xi, 7, = — — — (3.13) 

Xi+i — Aj 

From now on we will always assume <5j > 0, as we demand that X(x) = ^iS) 1 -^Q%( a? ) t* e a 
homeomorphism. We also have 

det ^ = -5 i _ 1 rf i (7 i _ 1 - 7i ) 2 . (3.14) 
y 

Proposition 3.4. T/ie mass matrix Mpj(y,X) is non-singular almost everywhere (as a function of 
the yi 's and Xi 's) and singular iff Ji-i = Ji for some i. 

Proof. We are going to compute the determinant of Mn{v,X) by transforming (13. IIP into a block 
upper triangular form by zeroing the blocks Bi below the diagonal. Let us start with the block B\. 
We use linear combinations of the first two rows of the mass matrix to zero the elements of the 
block B\ below the diagonal. Suppose 70 =71. Then it is easy to see that the first two rows of the 
mass matrix are not linearly independent, so the determinant of the mass matrix is zero. Assume 
7o 7^ 7i ■ Then by f|3. 14j) the block A\ is invertible. We multiply the first two rows of the mass 
matrix by B\A^ 1 and subtract the result from the third and fourth rows. This zeroes the block B\ 
below the diagonal and replaces the block A2 by 

C 2 = A 2 - B 1 A^ 1 B 1 . (3.15) 
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We now zero the block B2 below the diagonal in a similar fashion. After re — 1 steps of this procedure 
the mass matrix is transformed into 



/ Ci B 1 

C2 Bi 



Cn B r , 



B n A 



n+1 



V 



Bn-i 
Bn-i An J 



(3.16) 



In a moment we are going to see that C n is singular iff 7 n _i = j n and in that case the two rows 
of the matrix above that contain C n and B n are linearly dependent, thus making the mass matrix 
singular. Suppose "y n -i / 7™, so that C n is invertible. In the next step of our procedure the block 
A n+ \ is replaced by 

C n+ i = A n+ i — B n C~ l B n . (3.17) 
Together with the condition C\ = A\ this gives us a recurrence. By induction on n we find that 



C n 



\5n-l + — \&n-lln-l ~ \&nln 



(3.18) 



and 



detCi = — (5j_i(5i(7i_i - ji) 2 , 



(3.19) 



which justifies our assumptions on the invertibility of the blocks C{. We can now express the 
determinant of the mass matrix as detCi • ... • detC^r. The final formula is 



det M N (y,X) = °^ ^T~T^ (7o ~ 7i) 2 ---(7jV-i ~7n) 2 - 



(3.20) 



We see that the mass matrix becomes singular iff = ji for some i and this condition defines a 
measure zero subset of M. 2N . 

□ 



Remark I. This result shows that the finite-dimensional system described by the semi-discrete 
Lagrangian (|3,8p is non-degenerate almost everywhere. This means that, unlike in the continuous 
case, the Euler-Lagrange equations corresponding to the variations of the y^'s and Xi's are inde- 
pendent of each other (almost everywhere) and the equations corresponding to the X^s are in fact 
necessary for the correct description of the dynamics. This can also be seen in a more general 
way. Owing to the fact we are considering a finite element approximation, the semi-discrete action 
functional Sn is simply a restriction of S, and therefore formulas (13. ip still hold. The corresponding 
Euler-Lagrange equations take the form 
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x 



X 



Figure 3.1: Left: If 7&_i 7^ 7^, then any change to the middle point changes the local shape of 
(f)(X,t). Right: If ~fk-i = 7fc, then there are infinitely many possible positions for (X k ,y k ) that 
reproduce the local linear shape of cf>(X,t). 



5iS[ip,X] -6<p(x,t) = 0, (3.21) 
6 2 S[<p,X] -SX(x,t) = 0, 

which must hold for all variations 8tp(x, t) = Yl^=i 3yi{t) r li{ x ) an d 6X(x, t) = X^iLi SXi(t)rji(x). Since 
we are working in a finite dimensional subspace, the second equation now does not follow from the 
first equation. To see this, consider a particular variation 5X(x,t) = 8X k (t)rj k (x) for some k, where 
5Xk 7k 0. Then we have 



-jk-l SX k (t) m(x), if Xk-i <x<x k , 
~lk SX k (t) r) k (x), ifx k <x< x k+1 , (3.22) 
0, otherwise, 



which is discontinuous at x = x k and cannot be expressed as Y2i=i $yi(i) r ]i( x ) f° r an y ^Ui{t)i unless 
7fc— l = Ik- Therefore, we cannot invoke the first equation to show that 52S[np,X] ■ 5X(x,t) = 0. 
The second equation becomes independent. 

Remark II. It is also instructive to realize what exactly happens when j k -i = 7fe- This means 
that locally in the interval [X k _i, X k+ i] the field (p(X,t) is a straight line with the slope 7^. It also 
means that there are infinitely many values (X k ,y k ) that reproduce the same local shape of <f)(X, t). 
This reflects the arbitrariness of X(x,t) in the infinite-dimensional setting. In the finite element 
setting, however, this holds only when the points (X k -i, y k -i), (X k ,y k ) and (X k+ i,y k+ i) line up. 
Otherwise any change to the middle point changes the shape of <p(X,t). See Figure l3~Tl 

3.4 Existence and uniqueness of solutions 

Since the Legendre Transform (I3.10p becomes singular at some points, this raises a question about 
the existence and uniqueness of the solutions to the Euler-Lagrange equations (|3.2ip . In this section 
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we provide a partial answer to this problem. We will begin by computing the Lagrangian symplectic 
form 



N 



Q N = ~^2 dyi A dpi + dXj A dSi, 



(3.23) 



i=i 



where pi and Si are given by ([3.9p . For notational convenience we will collectively denote q = 
(yi, X u ...,y N ,X N ) T and q = ( yi , X 1: y N , X N ) T . Then in the ordered basis (^-, 
the symplectic form can be represented by the matrix 



A N (q,q) M N (q) 



-M N (q) 

where the 2N x 2iV block Apj(q, q) has the further block tridiagonal structure 



/ r x Ax 



An(q,> 



\ 



-Af r 2 a 2 
-Kl r 3 a 3 



Atv-i 



with the 2x2 blocks 



(3.24) 



(3.25) 



Vi+l—Vi-l Xj--i+2Xj 



7i-l + 



2Xi 4-Xi. 



-a 



A; 



3 ' 3 

2. 

Vi+i-yi _|_ 2Xi+x i+1 1 



7i-l 



7i 



6 1 3 
In this form, it is easy to see that 



6 3 '* 

2 'i 



(3.26) 



det^7v(g,g) = f detMjv(g) 

so the symplectic form is singular whenever the mass matrix is. 

The energy corresponding to the Lagrangian (j3.8H can be written as 



(3.27) 



1 r x k+i / 

En(q,Q) = -jQ r M N {q)q + 2_ J / R \lk, VkVk(x) + y k+1 rj k+1 ( 

k=o Xk 



.r||— — — dx. (3.25 
Ax 



In the chosen coordinates, gLEjv can be represented by the row vector dE^ = {dEpj / dq\, ...,dE^[/dq2N) 
It turns out that 
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where the vector £ has the following block structure 



(3.29) 



(3.30) 



Cn J 



Each of these blocks has the form = 1, Ck,2) T ■ Through basic algebraic manipulations and 
integration by parts, one finds that 

t _ ilk+i(2X k+ i + X k ) + y k (X k+1 - X k -\) - Vk-i{Xk + 2X fe _x) 
- g 

Xf + X k X k ^i + + x k+1 x k + xf 

+ o 7fc-i o 7fc 



and 



+ 



1 f Xk OR 



Ax Jx k _i d</>x 
1 r Xk +^ dR 
Ax _L, 84>x 



'X k 



+ 



ik,ykVk(x) + yk+iVk+i(x)^j dx 
l f Xk 



lk-l L 
1 



R(ik-i,yk) - J R{ik-i,yk-ir]k-i{x) + ykVk(x^j dx 



l f Xk+i / \ 

R {lk,yk)-j^ j R(j k ,ykVk{x) + y k+ ir] k+1 (x)) dx 



Xk 



(3.31) 



ik,2 — 



ik-i + iik-iik - Mk+i - y k+ i 



X\ + X k X k _i + x\_ x 2 + X k+ iX k + x| 



Ik 



7fc -i f Xk &R ( \ 
^7 7 ^l^- 1 '^- 1 *- 1 ^) + ykVk{x)J dx 

Ik I 
Ax Jx k 



(■Xk+1 QH , s 

+ -t— / 77— (jfc, yfe%(») + yk+iVk+i(x)j dx 



+ 



I fXk 



, , Rnk-uyk-iVk-iix) + ykVk(x)) dx 

^ x Jx k -i 

~~ ~Kx / R {^k,ykVk(x) + y k +m+i{x)) dx. 
We are now ready to consider the generalized Hamiltonian equation 



iz^N = dEjy, 
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(3.32) 



(3.33) 



which we solve for the vector field Z = X^i=i a i d/dqi + /3j d/dqi. In the matrix representation this 
equation takes the form 

nl(q,q)-( a A=dE%(q,q). (3.34) 

Equations of this form are called (quasilinear) implicit ODEs (see [44], |46|). If the symplectic form 
is nonsingular in a neighborhood of (q(°\ q^), then the equation can be solved directly via 

Z=[nl(q,q)]- 1 dEl(q,q) 

to obtain the standard explicit ODE form and standard existence/uniqueness theorems (Picard's, 
Peano's, etc.) of ODE theory can be invoked to show local existence and uniqueness of the flow of 
Z in a neighborhood of (q^°\q^). If, however, the symplectic form is sing ular at (q^,^), then 
there are two possibilities. The first case is 

dE%(q(°\qW) Range ^(0)^(0)) (3 M ) 

and it means there is no solution for Z at (q^°\q^)- This type of singularity is called an algebraic 
one and it leads to so called impasse points (see |39|-|44]. [46j). 

The other case is 

dEW°\Q {0) ) 6 Range n T N (q^,q^) (3.36) 

and it means that there exists a nonunique solution Z at (q^°\ q^)- This type of singularity is called 
a geometric one. If (q(°\q^) is a limit of regular points of (I3.34P (i.e. points where the symplectic 
form is nonsingular), then there might exist an integral curve of Z passing through (q^°\ q^)- See 
[39j-[44j, [46j for more details. 

Proposition 3.5. The singularities of the symplectic form Qi\[(q,q) are geometric. 

Proof Suppose that the mass matrix (and thus the symplectic form) is sing ular at (<? (0) , g (0) ). Using 
the block structures (I3.24P and (13.29|) we can write (13.34P as the system 

-A N (q^,q^)a-M N (q^)f3 = ^ 

M N (qW)a = M N {qW)q ( -°\ (3.37) 

The second equation implies that there exists a solution a = q(°\ In fact this is the only solution 
we are interested in, since it satisfies the second order condition: the Euler-Lagrange equations 
underlying the variationl principle are second order, so we are only interested in solutions of the 
form Z = X)i=i Qi 9/dqi + Pi d/dqi. The first equation can be rewritten as 

M N (q^)P = -t-h N (q iQ \q^)q^. (3.38) 
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Since the mass matrix is singular, we must have 7fe_i = jk f° r some k. As we saw in Section 13,31 
this means that the two rows of the k th 'block row' of the mass matrix (i.e., the rows containing 
the blocks B^-i, Ak and B^) are not linearly independent. In fact we have 

(£^-1)2* = -7fc(-Bfc-i)i*, {Ak)2* = -7k(Ak)i*, {Bk)2* = -7k(Bk)u, (3.39) 

where a m * denotes the m th row of the matrix a. Equation (|3,38p will have a solution for (3 iff 
the RHS satisfies a similar scaling condition in the the k th 'block element'. Using formulas (|3.26p . 
(|3.3ip and ([3.32p . we show that — £ — An q^ indeed has this property. Hence, dE N (q^°\q^) 6 
Range f2&(gW,?W) and (g (0) ,g (0 )) is a geometric singularity. Moreover, since jk-i = 7fe defines a 
hypersurface in R 2N x R 2N , (<? (0) ,g (0) ) is a limit of regular points. 

□ 

Remark I. Numerical time integration of the semi-discrete equations of motion (13. 34ft has to 
deal with the singularity points of the symplectic form. While there are some numerical algorithms 
allowing one to get past singular hypersurfaces (see |44j), it might not be very practical from the 
application point of view. Note that, unlike in the continuous case, the time evolution of the 
meshpoints Xj's is governed by the equations of motion, so the user does not have any influence 
on how the mesh is adapted. More importantly, there is no built-in mechanism that would prevent 
mesh tangling. Some preliminary numerical experiments show that the mesh points eventually 
collapse when started with nonzero initial velocities. 

Remark II. The singularities of the mass matrix (13. lip bear some similarities to the singularities 
of the mass matrices encountered in the Moving Finite Element method. In |35j and [36] the authors 
proposed introducing a small 'internodal' viscosity which penalizes the method for relative motion 
between the nodes and thus regularizes the mass matrix. A similar idea could be applied in our case: 
one could add some small e kinetic terms to the Lagrangian (|3.8p in order to regularize the Legendre 
Transform. In light of the remark made above, we did not follow this idea further and decided to 
take a different route instead, as described in the following sections. However, investigating further 
similarities between our variational approach and the Moving Finite Element method might be 
worthwhile. There also might be some connection to the r-adaptive method presented in |51j : the 
evolution of the mesh in that method is also set by the equations of motion, although the authors 
considered a different variational principle and different theoretical reasoning to justify the validity 
of their approach. 

3.5 Constraints and adaptation strategy 

As we saw in Section T3.41 upon discretization we lose the arbitrariness of X(x,t) and the evolution 
of Xi(t) is governed by the equations of motion, while we still want to be able to select a desired 
mesh adaptation strategy, like (|2.25p . This could be done by augmenting the Lagrangian (|3.8p 
with Lagrange multipliers corresponding to each constraint g{. However, it is not obvious that 
the dynamics of the constrained system as defined would reflect in any way the behavior of the 
approximated system (13 . 6|) . We are going to show that the constraints can be added via Lagrange 
multipliers already at the continuous level (|3.6p and the continuous system as defined can be then 
discretized to arrive at ([3.8p with the desired adaptation constraints. 
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3.5.1 Global constraint 

As mentioned before, eventually we would like to impose the constraints 



g i (y 1 ,...,y N ,X 1 ,...,X N ) = i = l,...,N (3.40) 

on the semi-discrete system (|3.8p . Let us assume that g : M. 2N — > M. N , g = (gi, gN) T is C 1 
and is a regular value of g, so that (|3.40p defines a submanifold. To see how these constraints 
can be introduced at the continuous level, let us select uniformly distributed points x% = i • Ax, 
i = 0, N + 1, Ax = X max /(N + 1) and demand that the constraints 

g i (<p(x 1 ,t),...,ip(x N ,t),X(x 1 ,t),...,X(x N ,t)\=0, i = l,...,N (3.41) 

be satisfied by (p(x,t) and X(x,t). One way of imposing these constraints is solving the system 



6xS[(p, X] ■ 5<p(x, t)=0 for all 8<p{x, t), (3.42) 
gi(p(xi,t),...,(p(x N ,t),X(xi,t),...,X(x N ,t)j = 0, i= l,...,N. 

This system consists of one Euler-Lagrange equation that corresponds to extremizing S with respect 
to ip (we saw in Section 13.11 that the other Euler-Lagrange equation is not independent) and a set 
of constraints enforced at some pre-selected points X{. Note, that upon finite element discretization 
on a mesh coinciding with the pre-selected points this system reduces to the approach presented 
in Section [2) we minimize the discrete action with respect to the y^s only and supplement the 
resulting equations with the constraints (|3.40D . 

Another way that we want to explore consists in using Lagrange multipliers. Define the auxiliary 
action functional 



N »T 

^ ^ / J- max / \ 

Sc[<f, X, Afc] = S[ip, X] - / \i(t)-gi\tp(xi,t),...,tp(xN,t),X(xi,t),...,X(xN,t)) dt. (3.43) 

i=i ^° 

We are going to assume that the Lagrange multipliers Aj(i) are at least continuous in time. According 
to the method of Lagrange multipliers, we seek the stationary points of Sc- This leads to the 
following system of equations 



~ ^ _^ r J-max dQ 

5\S\ip, X] ■ 6ip(x, t) — N N / Aj(i) —-^ 6ip(xj, t) dt = for all 5tp(x, t), 

N N ,T rs 

f J- max fin ■ 

S 2 S[(p,X] -SX(x,t)- y^Yl / X i (t)-^5X{x j ,t)dt = for all SX(x,t), 
i=i j=i Jo dX i 

g i (<p(x 1 ,t),...,(p{x N ,t),X(x 1 ,t),...,X(x N ,t)) =0, i = l,...,N, (3.44) 



where for clarity we suppressed writing the arguments of ^ and Jj^. 
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Equation (|3,42p is more intuitive, because we directly use the arbitrariness of X(x, t) and simply 
restrict it further by imposing constraints. It is not immediately obvious how solutions of (|3.42p 
and (I3.44p relate to each other. We would like both systems to be 'equivalent' in some sense, or at 
least their solution sets to overlap. Let us investigate this issue in more detail. 

Suppose (<p,X) satisfy (|3.42p . Then it is quite trivial to see that {ip,X, Ai, Ajv) such that 
Afc = satisfy (|3.44j) : the second equation is implied by the first one and the other equations 
coincide with those of (13.42p . At this point it should be obvious that system (I3.44P may have more 
solutions for <p and X than system (I3.42p . 

Proposition 3.6. The only solutions (</?, X, X%, Ajv) to (|3.44p that satisfy (|3.42p as well are those 
with Afc = for all k. 

Proof. Suppose (ip, X, Ai, Ajv) satisfy both (I3.42p and (|3.44j) . System (|3,42p implies that 5iS-5p = 
and 62 S ■ 5X = 0. Using this in system (|3.44|) gives 



N „T N 



/ dt 5<p(xj,t) Xj(t) = for all 5<p (x, t), 

j 1 i=i ,>!J > 

N »r 

I J- max f~)n 

V/ dt6X( Xj ,t) VAi(t) -7^- = for all SX(x,t). (3.45) 

T"W0 ~~1 uXj 



j=l" u i=l ' 3 



In particular, this has to hold for variations Sep and SX such that 5<p(xj,t) = 5X(xj,t) = v{t) ■ 5kj, 
where v{t) is an arbitrary continuous function of time. If we further assume that for all x E 
[0, X max ] the functions ip(x, .) and X(x, .) are continuous, both X/i=i ^i(^) an d X^i=i ^i(^) 
are continuous and we get 

Z? 5 ^(xi,t),...,<^(a; 7V ,t),X(xi,t),...,X( a ; i v,t)) • A(t) = (3.46) 



for all t, where A = (Ai, Ajv) and the iV x 2N matrix Do = jJt- is the derivative 

of 5. Since we assumed that is a regular value of g and the constraint g = is satisfied by tp and 
X, we have that for all t the matrix L>g has full rank — that is, there exists a nonsingular N x N 
submatrix S. Then the equation H T A(i) = implies A = 0. 

□ 

We see that considering Lagrange multipliers in (I3.43P makes sense at the continuous level. We 
can now perform a finite element discretization. The auxiliary Lagrangian Lc :QxGxWxZx 
W N — > R corresponding to (13.43P can be written as 



N 

Lc[<p,X,<pt,X t ,\k] = L[p,X,p t ,X t ] - y^Aj ■ gi(tp(xi),...,tp(x N ),X(xi),...,X(x N )\ (3.47) 

i=i 

where L is the Lagrangian of the unconstrained theory and has been defined by (13 . 2|) . Let us 
choose a uniform mesh coinciding with the pre-selected points X{. As in Section [3.21 we consider 
the restriction L CN = L \q nX g n xW n xZ n x^ and we g et 
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N 

LcN(yi, Xj,y k , X h A m ) = L N (yi,Xj,y k ,Xi) - ^ Aj ■ g^, ...,y N , X 1 , X N ). (3.48) 

i=l 

We see that the semi-discrete Lagrangian Lqn is obtained from the semi-discrete Lagrangian Ljv by 
adding the constraints g^ directly at the semi-discrete level, which is exactly what we set out to do at 
the beginning of this section. However, in the semi-discrete setting we cannot expect the Lagrange 
multipliers to vanish for solutions of interest. This is because there is no semi-discrete counterpart 
of Proposition 13.61 On one hand, the semi-discrete version of (|3.42p (that is, the approach presented 
in Section[2]) does not imply that b^S-bX = 0, so the above proof will not work. On the other hand, 
if we supplement (|3.42p with the equation corresponding to variations of X, then the finite element 
discretization will not have solutions, unless the constraint functions are integrals of motion of the 
system described by L]\f(yi,Xj,yi t ,Xi), which generally is not the case. Nonetheless, it is reasonable 
to expect that if the continuous system (|3.42p has a solution, then the Lagrange multipliers of the 
semi-discrete system (|3.48p should remain small. 

Defining constraints by Equations ([3.4ip allowed us to use the same finite element discretization 
for both L and the constraints, and to prove some correspondence between the solutions of (I3.42P 
and (|3,44p . However, constraints (|3,4ip are global in the sense that they depend on the values of 
the fields ip and X at different points in space. Moreover, these constraints do not determine unique 
solutions to (|3.42p and f|3.44j) , which is a little cumbersome when discussing multisymplecticity (see 
Section |4]). 

3.5.2 Local constraint 

In Section \2. 41 we discussed how some adaptation constraints of interest can be derived from certain 
partial differential equations based on the equidistribution principle, for instance equation (|2.24p . 
We can view these PDEs as local constraints that only depend on pointwise values of the fields 
<p, X and their spatial derivatives. Let G = G((p, X, (p x , X x , (f xx , X xx , ...) represent such a local 
constraint. Then, similarly to f|3.42|) . we can write our control-theoretic strategy from Section [2] as 

<5i5[<^ ; X] • 5<p(x, = for all Sip(x, t), (3.49) 
G(cp,X, (p x ,X x x xx , ...) — 0. 

Note that higher order derivatives of the fields may require the use of higher degree basis functions 
than the ones in ([2.6p . or of finite differences instead. 

The Lagrange multiplier approach consists in defining the auxiliary Lagrangian 

£>c[<P, X, <pt, X t , A] = L[ip, X, ipt,X t ] — / \(x) ■ G(cp,X,ip x ,X x ,(p xx ,X xx , ...)dx. (3.50) 

J o 

Suppose that the pair (ip,X) satisfies (|3,49p . Then, much like in Section [3.5. 1\ one can easily check 
that the triple (ip, X, A = 0) satisfies the Euler-Lagrange equations associated with (13.50|) . However, 
an analog of Proposition 13.61 does not seem to be very interesting in this case, therefore we are not 
proving it here. 
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Introducing the constraints this way is convenient, because the Lagrangian (|3,50p then represents 
a constrained multisymplectic field theory with a local constraint, which makes the analysis of 
multisymplecticity easier (see Section 2]). The disadvantage is that discretization of (I3.50P requires 
mixed methods. We are going to use the linear finite elements (|2.6p to discretize L[<p,X,<pt,Xt], 
but the constraint term will be approximated via finite differences. This way we again obtain the 
semi-discrete Lagrangian (|3.48p . where gi represents the discretization of G at the point x = x%. 

In summary, the methods presented in Section 13.5.11 and Section 13.5.21 both lead to the same 
semi-discrete Lagrangian, but have different theoretical advantages. 

3.6 DAE formulation of the equations of motion 

The Lagrangian (|3.48|) can be written as 

L CN (q, q, A) = \q T M N {q) q - R N (q) - X T g(q), (3.51) 

where 



R N{q) = ^2f R(ik, Vkm(x) + yk+iVk+i(xf) Xh+l ^ Xk dx - ( 3 - 52 ) 



The Euler-Lagrange equations thus take the form 



u, 



M N (q)u = f(q,u)-Dg(qy A, 

g(q) = 0, (3.53) 



where 



2N 



8R N x ^ (I djM^ij 8{M N ) kt x 



dqk .^H, ^2 dq k dq 

System f|3.53j) is to be solved for the unknown functions q(t), u(t) and \(t). This is a DAE system 
of index 3, since we are lacking a differential equation for X(t) and the constraint equation has to 
be differentiated three times in order to express A as a function of q, u and A, provided that certain 
regularity conditions are satisfied. Let us determine these conditions. Differentiate the constraint 
equation with respect to time twice to obtain the acceleration level constraint 

Dg(q)u = h(q,u), (3.55) 

where 

2N n2 



3 
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We can then write (|3,55p and the second equation of (|3.53p together as 



M N (q) Dg(q) T \( u \ _ ( f(q,u) \ 

Dg(q) ){\ )-{h(q,u) J' ^ 

If we could solve this equation for u and A in terms of q and u, then we could simply differentiate 
the expression for A one more time to obtain the missing differential equation, thus showing system 
(|3.53p is of index 3. System (|3.57p is solvable if its matrix is invertible. Hence, for system (|3.53p to 
be of index 3 the following condition 

has to be satisfied for all q or at least in a neighborhood of the points satisfying g(q) = 0. Note 
that with suitably chosen constraints this condition allows the mass matrix to be singular. 

We would like to perform time integration of this mechanical system using the symplectic (vari- 
ational) Lobatto IIIA-IIIB quadratures for constrained systems (see |20j . [22], [26] . [27], [34] ). How- 
ever, due to the singularity of the Runge-Kutta coefficient matrices (a«) and iflij) for the Lobatto 
IIIA and IIIB schemes, the assumption (I3.58P does not guarantee that these quadratures define 
a unique numerical solution: the mass matrix would need to be invertible. To circumvent this 
numerical obstacle we resort to a trick described in [27] . We embed our mechanical system in a 
higher dimensional configuration space by adding slack degrees of freedom r and r and form the 
augmented Lagrangian Lj^ by modifying the kinetic term of L^r to read 

£*to,r,4,o -\(e r-)-( ° 3 f )■(*)" «»<*>■ < 3 - 59 > 

Assuming (I3.58p . the augmented system has a non-singular mass matrix. If we multiply out the 
terms we obtain simply 

L$(q, r, q, r) = L N (q, q) + r T Dg(q) q. (3.60) 

This formula in fact holds for general Lagrangians, not only for (13. 8D . In addition to g(q) = we 
further impose the constraint r = 0. Then the augmented constrained Lagrangian takes the form 

Lcn(Q> r > 9i f, A, fj) = L N (q, q) + r T Dg(q) q - \ T g{q) - [i T r. (3.61) 
The corresponding Euler-Lagrange equations are 



M N (q)u + Dg(qfw 
Dg{q) u 

9(q) 



w, 

f(q,u) 

h(q,u) 

0, 

0. 



Dg(q) T X, 
Mi 



(3.62) 
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It is straightforward to verify that r(t) = 0, w(t) = 0, fi(t) = is the exact solution and the 
remaining equations reduce to (|3,53p . that is, the evolution of the augmented system coincides with 
the evolution of the original system, by construction. The advantage is that the augmented system 
is now regular and we can readily apply the Lobatto IIIA-IIIB method for constrained systems 
to compute a numerical solution. It should be intuitively clear that this numerical solution will 
approximate the solution of (|3,53p as well. What is not immediately obvious is whether a variational 
integrator based on (13.60|) can be interpreted as a variational integrator based on Lj\r. This can 
be elegantly justified with the help of exact constrained discrete Lagrangians. Let TV C Qn x Gn 
be the constraint submanifold defined by g(q) = 0. The exact constrained discrete Lagrangian 
L% E : TV x TV — ► K is defined by 

/■At 

L% E {qW,qW) = / L N {q(t),q(t))dt, (3.63) 
J o 

where q(t) is the solution to the constrained Euler-Lagrange equations (|3.53|) such that it satisfies 
the boundary conditions q(0) = <? (1) and q(At) = q^ . Note that TV x {0} C (Qn x G n ) x R n is 
the constraint submanifold defined by g(q) = and r = 0. Since necessarily rW = A 2 ^ = 0, we can 
define the exact augmented constrained discrete Lagrangian L N ' : : TV x TV — > R by 

r-At 

L% C > E (q< 1 \q®) = / Li(q(t),r(t),q(t),r(t))dt, (3.64) 
J o 

where q(t), r(t) are the solutions to the augmented constrained Euler-Lagrange equations (13. 62^ 
such that the boundary conditions q(Q) = q^\ q(At) = q^ and r(0) = r(At) = are satisfied. 

Proposition 3.7. The exact discrete Lagrangians Lfj ,C ' E and L N ' E are equal. 

Proof. Let q(t) and r(t) be the solutions to fl3.62[> such that the boundary conditions q(0) = q^\ 
q(At) = qW and r(0) = r(At) = are satisfied. As argued before, we in fact have r(t) = and q(t) 
satisfies (I3.53P as well. By (I3.60P we have 

L N (q(t),r(t),q(t),f(t)) = L N (q{t),q(t)) 
for all t G [0, At], and consequently L N ' C,E = L N ' E . 

□ 

This means that any discrete Lagrangian : (Qn x Gn) x R n x (Qn x Gn) x R n — > R that 
approximates L N ' C ' E to order s also approximates L N ' E to the same order, that is, a variational 
integrator for (13.62H . in particular our Lobatto IIIA-IIIB scheme, is also a variational integrator for 
(13331) . 

Backward error analysis. The advantage of the Lagrange multiplier approach is the fact that 
upon spatial discretization we deal with a constrained mechanical system. Backward error analysis 
of symplectic / variational numerical schemes for such systems shows that the modified equations also 
describe a constrained mechanical system for a nearby Hamiltonian (see |20|). Therefore, we expect 
the Lagrange multiplier strategy to demonstrate better performance in terms of energy conservation 
than the control-theoretic strategy. The Lagrange multiplier approach makes better use of the 
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geometry underlying the field theory we consider, the key idea being to treat the reparametrization 
field X(x,t) as an additional dynamical degree of freedom on equal footing with tp(x,t). 

4 Multisymplectic field theory formalism 

In Section [2] and Section [3] we took the view of infinite dimensional manifolds of fields as configuration 
spaces and presented a way to construct space-adaptive variational integrators in that formalism. 
We essentially applied symplectic integrators to semi-discretized Lagrangian field theories. In this 
section we show how r-adaptive integrators can be described in the more general framework of 
multisymplectic geometry. In particular we show that some of the integrators obtained in the 
previous sections can be interpreted as multisymplectic variational integrators. Multisymplectic 
geometry provides a covariant formalism for the study of field theories in which time and space are 
treated on equal footing, as a conseqence of which multisymplectic variational integrators allow for 
more general discretizations of spacetime, such that, for instance, each element of space may be 
integrated with a different timestep (see [30J). For the convenience of the reader, below we briefly 
review some background material and provide relevant references for further details. We then 
proceed to reformulate our adaptation strategies in the language of multisymplectic field theory. 

4.1 Background material 

Lagrangian mechanics and Veselov-type discretizations 

Let Q be the configuration manifold of a certain mechanical system and TQ its tangent bundle. 
Denote the coordinates on Q by q l , and on TQ by (q l ,q l ), where i = 1,2, ...,n. The system 
is described by defining the Lagrangian L : TQ — > R and the corresponding action functional 
(fjf^i)] = f^L(q l (t),q l (t)) dt. The dynamics is obtained through Hamilton's principle, which seeks 
the curves q(t) for which the functional S[q(i)] is stationary under variations of q(t) with fixed 
endpoints, i.e. we seek q{t) such that 



dS[q(t)] ■ 5q(t) = - 



S[q t (t)] = 



(4.1) 



for all 8q(t) with 6q(a) = Sq(b) = 0, where q e (t) is a smooth family of curves satisfying q$ 
^\ e —oQe = Sq. By using integration by parts, the Euler-Lagrange equations follow as 



and 



dq i dtdq 1 °" ^'^ 

The canonical symplectic form Q on T*Q, the 2n-dimensional cotangent bundle of Q, is given by 
0, = dq l Adpi, where summation over i is implied and (q l ,Pi) are the canonical coordinates on T*Q. 
The Lagrangian defines the Legendre transformation FL : TQ — > T*Q, which in coordinates is 
given by (q l ,Pi) = (q l , |pr)- W e then define the Lagrange 2-form on TQ by pulling back the canonical 
symplectic form, i.e. VIl = FL*0. If the Legendre transformation is a local diffeomorphism, then 
CIl i s a symplectic form. The Lagrange vector field is a vector field on TQ that satisfies 
XejQl = dE, where the energy E is defined by E(v q ) = ¥L(v q ) • v q — L(v q ) and j denotes the 
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interior product, i.e. the contraction of a differential form with a vector field. It can be shown that 
the flow Ft of this vector field preserves the symplectic form, that is, F^VLl = £1^. The flow Ft is 
obtained by solving the Euler-Lagrange equations (14. 21 . 

For a Veselov-type discretization we essentially replace TQ with Q x Q, which serves as a 
discrete approximation of the tangent bundle. We define a discrete Lagrangian L d as a smooth map 
Ld '■ Q x Q — > R and the corresponding discrete action S = Ylk=o Ld(q k ,q k +i). The variational 
principle now seeks a sequence qo, qi, qN that extremizes S for variations holding the endpoints 
qo and qjy fixed. The Discrete Euler-Lagrange equations follow 

D 2 L d (q k ^ 1 ,q k ) + D 1 L d (q k , q k+l ) = 0. (4.3) 

This implicitly defines a discrete flow F : Q x Q — > Q x Q such that F(q k ^i,q k ) = (q k , q k +i)- One 
can define the discrete Lagrange 2-form on Q X Q by ojl = q dq}) A dq^, where (qQ,q{) denotes 
the coordinates on Q x Q. It then follows that the discrete flow F is symplectic, i.e. F*ujl = ^L- 

Given a continuous Lagrangian system with L : TQ — > M one chooses a corresponding discrete 
Lagrangian as an approximation L^qk^qk+i) ~ ft* +1 L(q(t), q(t)} dt, where q(t) is the solution 
of the Euler-Lagrange equations corresponding to L with the boundary values q{t k ) = q k and 
Q^fc+l) = Qk+l- 

For more details regarding Lagrangian mechanics, variational principles, and symplectic geom- 
etry, see [33J. Discrete Mechanics and variational integrators are discussed in [34j. 

Multisymplectic geometry and Lagrangian field theory 

Let X be an oriented manifold representing the (n + l)-dimensional spacetime with local coordi- 
nates x , . . . , x n ) = (t,x), where x° = t is time and (x 1 ,...,x n ) = x are space coordinates. 
Physical fields are sections of a configuration fiber bundle ttxy '■ Y — > X, that is, continuous maps 
(f) : X — > Y such that ttxy 4> = id^. This means that for every (t, x) G X, <p(t,x) is in the 
fiber over (t,x), which is Yu x ^ = ir XY ((t,x)). The evolution of the field takes place on the first 
jet bundle J 1 !^, which is the analog of TQ for mechanical systems. J 1 !^ is defined as the affine 
bundle over Y such that for y G Y(t,x) the fiber JyY consists of linear maps 1? : Tr tjX )X — > T y Y 
satisfying the condition Tttxy $ = '^Tu^x- The local coordinates (x^,y a ) on Y induce the co- 
ordinates (x^, y a , v a ^) on J l Y . Intuitively, the first jet bundle consists of the configuration bundle 
Y, and of the first partial derivatives of the field variables with respect to the independent vari- 
ables. Let <p(x°, . . . , x n ) = (ar, . . . , x n , y , . . . , y m ) in coordinates and let u a = y a ^ = dy a /dx fi 
denote the partial derivatives. We can think of J l Y as a fiber bundle over X. Given a sec- 
tion cj) '■ X — > Y, we can define its first jet prolongation j 1 ^ : X — > J l Y , in coordinates 
given by j 1 (j)(x°,x 1 , . . . , x n ) = x 1 , . . . , x n , y , . . . , y m , y l Q , . . . , y m n ), which is a section of the 
fiber bundle J X Y over X. For higher order field theories we consider higher order jet bundles, 
defined iteratively by J 2 Y = J 1 (J 1 !") and so on. The local coordinates on J 2 Y are denoted 
y a , v a ^ : w a ^, Knv)- The second jet prolongation j 2 cf) : X — > J 2 Y is given in coordinates by 
j 2 Hx») = (x»,y a ,y\,y%,y\ ;1/ ). 

Lagrangian density for first order field theories is defined as a map C : J l Y — > M. The 
corresponding action functional is S[4>] = J u C^ 1 4>) d n+1 x , where U C X. Hamilton's principle 
seeks fields cj>(t, x) that extremize S, that is 
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S[4 O 0] = (4.4) 

dA A=0 

for all t]y that keep the boundary conditions on dlA fixed, where riy '■ Y — > Y is the flow of a 
vertical vector field V on Y. This leads to the Euler-Lagrange equations 

^"■sfl^T 11 (4 - 5) 

Given the Lagrangian density C one can define the Cartan (n + l)-form 0£ on J 1 !^, in local 
coordinates given by 0£ = ^-dy a A + (£ — ■^j=-v a fl )d n+1 x, where d n x fl = 9 M The 

multisymplectic (n + 2)-form is then defined by £lc = —d@c- Let V be the set of solutions of 
the Euler-Lagrange equations, that is, the set of sections <j) satisfying (|4,4p or (|4.5p . For a given 
<f> G V, let T be the set of first variations, that is, the set of vector fields V on J^Y such that 
(t, x) — > r]y o (j)(t, x) is also a solution, where r/y is the flow of V . The multisymplectic form formula 
states that if 6 G V then for all V and PF in J-, 



f (j 1 <f>)*{j 1 V^j 1 Wjn c )=0, (4.6) 

JdU 

where j l V is the jet prolongation of V, that is, the vector field on J^Y in local coordinates given 
by fV = (V,V a , + ^rv\ - v a v ^-), where V = (V^,V a ) in local coordinates. The mul- 
tisymplectic form formula is the multisymplectic counterpart of the fact that in finite-dimensional 
mechanics, the flow of a mechanical system consists of symplectic maps. 

For a fc th -order Lagrangian field theory with the Lagrangian density C : J k Y — > R, analo gous 
geometric structures are defined on J 2k ~ 1 Y. In particular, for a second-order field theory the 
multisymplectic (n + 2)-form Qc is defined on J 3 Y and a similar multisymplectic form formula can 
be proven. If the Lagrangian density does not depend on the second order time derivatives of the 
field, it is convenient to define the subbundle JqY C J 2 Y such that JqY = {■& £ J 2 Y \ Kq = 0}. 

For more information about the geometry of jet bundles, see |47| . The multisymplectic formalism 
in field theory is discussed in [19J. The multisymplectic form formula for first-order field theories is 
derived in |31| . and generalized for second-order field theories in |28j. Higher order field theory is 
considered in [18]. 



Multisymplectic variational integrators 

Veselov-type discretization can be generalized to multisymplectic field theory. We take X = Z x Z = 
{(j, i)}, where for simplicity we consider dimAf = 2, i.e. n = 1. The configuration fiber bundle is 
Y = X x & for some smooth manifold & . The fiber over ( j, i) G X is denoted Yji and its elements 
Uji. A rectangle □ of X is an ordered 4-tuple of the form □ = ((J, i), (j, (j + 1, + i)) = 

(D 1 , D 2 , D 3 , D 4 ). The set of all rectangles in X is denoted X u . A point (j, i) is touched by a rectangle 
if it is a vertex of that rectangle. Let IA C X . Then (J, i) G IA is an interior point of U if £Y contains 
all four rectangles that touch The interior intW is the set of all interior points of U. The 

closure clU is the union of all rectangles touching interior points of U. The boundary of U is defined 
by dU = (U n clW)\int W. A section of Y is a map 4> : U C X — > Y such that i) G Yj-j. We can 
now define the discrete first jet bundle of Y as 
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J rl * r = {(yji,yji+i,yj+ii+i,yj+n) \ e ^, yji,yji+i,yj+n+i,yj+n e ^} 

= Af n x J^ 4 . (4.7) 

Intuitively, the discrete first jet bundle is the set of all rectangles together with four values assigned to 
their vertices. Those four values are enough to approximate the first derivatives of a smooth section 
with respect to time and space using, for instance, finite differences. The first jet prolongation of a 
section of Y is the map jV : X u -> J X Y defined by jV(P) = (□, 4>( al ), <AP 2 ), 0<P 3 ), 0<P 4 ))- 
For a vector field V on Y, let V^j be its restriction to Yji. Define a discrete Lagrangian L : J l Y — > R, 
L = ?/2 5 2/3, 2/4)) where for convenience we omit writing the base rectangle. The associated 

discrete action is given by 

s[<f>] = Y, L °3 1 <K n )- 

new 

The discrete variational principle seeks sections that extremize the discrete action, that is, mappings 
4>(j, i) such that 



d_ 

dX 



S[<f>x] = (4. 

A=0 



for all vector fields V on Y that keep the boundary conditions on dlA fixed, where (f)\(j,i) 
equations 



((j)(J,i)) and fY 1 is the flow of Vji on This is equivalent to the discrete Euler-Lagrange 



dL dL 

-K—{yju yji+i,yj+ii+i,yj+n) + -^—{yji-uyji, yj+ihyj+n-i)+ 

oyi oy 2 

dL dL . , . • 

+ a—{yj-ii-iiyj-u,yji,yji-i) + -^—{yj-ihyj-ii+uyji+uyji) = o ( 4 -9) 

for all £ mtU, where we adopt the convention <j)(j,i) = yj%. In analogy to the Veselov 

discretization of mechanics, we can define four 2-forms Q l L on J l Y , where I = 1,2,3,4 and Q} L + 
+ + 0^ = 0, that is, only three 2-forms of these forms are independent. The 4-tuple 
O?, Of,) is the discrete analog of the multisymplectic form Q^. We refer the reader to the 
literature for details, e.g. |31j . By analogy to the continuous case, let V be the set of solutions of 
the discrete Euler-Lagrange equations (|4.9p . For a given <p 6 V, let T be the set of first variations, 
that is, the set of vector fields V on ^Y defined similarly as in the continuous case. The discrete 
multisymplectic form formula then states that if 6 V then for all V and W in J-, 



E ( E [u^ru'vjj'wj^) 
□ v 1 



(□) = 0, (4.10) 



where the jet prolongations are defined to be 



j 1 V(y n i,y D 2,y D 3,y n 4) = (V D i(y D i), V D 2 (y n2 ), V D 3(y a s), Vb^yoO). (4.11) 
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The discrete form formula (|4,10p is in direct analogy to the multisymplectic form formula (|4,6p that 
holds in the continuous case. 

Given a continuous Lagrangian density C one chooses a corresponding discrete Lagrangian as 
an approximation L(y n i,y n 2,?/ n 3, y\jt) ~ ° j 4>dx dt, where □ is the rectangular region of the 
continuous spacetime that contains □ and (f)(t, x) is the solution of the Euler-Lagrange equations 
corresponding to C with the boundary values at the vertices of □ corresponding to y C 2, y C 3, 
and j/q4. 

The discrete second jet bundle J 2 Y can be defined by considering ordered 9-tuples 



ffl = {(j 1), (j - (J - 1,* + 1), 0>' - 1), 

(j,i),CM + i),(j + M- 1),0' + 1,0. (i + M + i)) 

= (ffl 1 , ffl 2 , ffl 3 , ffl 4 , ffl 5 , ffl 6 ,ffl 7 , ffl 8 , ffl 9 ) (4.12) 
instead of rectangles □, and the discrete subbundle JqY can be defined by considering 6-tuples 



m = - l), (j,0, + + + + 1,0, (j + M-1)) 
= (m 1 ,m 2 ,m 3 ,m 4 ,[i] 5 ,[i] 6 ). (4.13) 

Similar constructions then follow and a similar discrete multisymplectic form formula can be derived 
for a second order field theory. 

Multisymplectic variational integrators for first order field theories are introduced in [31j , and 
generalized for second-order field theories in [28J. 



4.2 Analysis of the control-theoretic approach 
Continuous setting 

We are now going to discuss a multisymplectic setting for the approach presented in Section [2j 
Let the computational spacetime be X = R x R with coordinates (t, x) and consider the trivial 
configuration bundle Y = X x R with coordinates (t,x,y). Let IA = [0, T max \ x [0,X max ] and 
let our scalar field be represented by a section tp : U — > Y with the coordinate representation 
<p(t,x) = (t,x,(p(t,x)). Let (t,x,y,Vt,v x ) denote local coordinates on J Y . In these coordinates 
the first jet prolongation of (p is represented by j tp(t,x) = (t,x,(p(t,x),ipt(t,x),(p x (t,x)). Then the 
Lagrangian density (12.3j) can be viewed as a mapping C : J^Y — > R. The corresponding action 
(|1.3p can now be expressed as 

S[0\= [ CtftydtAdx, (4.14) 

Just like in Section [21 let us for the moment assume that the function X : IA — > [0,X max ] is known, 
so that we can view C as being time and space dependent. The dynamics is obtained by extremizing 
S with respect to <p, that is, by solving for <p such that 

§[rfco0\ = O (4.15) 

A=0 



d 

dX 
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for all t]y that keep the boundary conditions on dlA fixed, where r/y '■ Y — > Y is the flow of a 
vertical vector field V on Y. Therefore, for an a priori known X(t, x) the multisymplectic form 
formula (14.61) is satisfied for solutions of (14. 15ft . 

Consider the additional bundle iixB ■ B = X X [0, X max ] — >• X whose sections X : U — > B 
represent our diffeomorphisms. Let X(t, x) = (t, x, X(t, x)) denote a local coordinate representation 
and assume X(t, .) is a diffeomorphism. Then define Y = Y © B. We have J k Y = J k Y © J k B. In 
Section 13.5.21 we argued that the moving mesh partial differential equation (12. 22j) can be interpreted 
as a local constraint on the fields <p, X and their spatial derivatives. This constraint can be repre- 
sented by a function G : J k Y — > R. Sections tip and X satisfy the constraint if G(j k <p,j k X) = 0. 
Therefore our control-theoretic strategy expressed in equations (|3.49|) can be rewritten as 



_d_ 

dX 



S[ V y O<p]=0, 

A=0 

G(j k <p,j k X) = 0, (4.16) 



for all rjy, similarly as above. Let us argue how to interpret the notion of multisymplecticity for this 
problem. Intuitively, multisymplecticity should be understood in a sense similar to Proposition 12.31 
We first solve the problem (14.16|) for tip and X , given some initial and boundary conditions. Then we 
substitute this X into the problem (|4.15p . Let V be the set of solutions to this problem. Naturally, 
tip G V . The multisymplectic form formula (I4.6P will be satisfied for all fields in V, but the constraint 
G = will be satisfied only for tip. 

Discretization 

Discretize the computational spacetime M x M by picking the discrete set of points tj = j ■ At, Xi = 
i • Ax, and define X = {(j, i) \j,i G Z}. Let X u and X m be the set of rectangles and 6-tuples in X, 
respectively. The discrete configuration bundle is Y = X x R and for convenience of notation let the 
elements of the fiber Yji be denoted by y\. Let U = {(j, i) \ j = 0, 1, . . . , M + 1, i = 0, 1, . . . , N + 1}, 
where Ax = X max /(N + 1) and At = T max /(M + 1). Suppose we have a discrete Lagrangian 
L : J l Y — > IR and the corresponding discrete action S that approximates (I4.14p . where we assume 
that X(t,x) is known and of the form (|2.7p . A variational integrator is obtained by solving 



d 

dX 



S[<px] = (4.17) 

A=0 



for a discrete section tip : U — > Y, as described in Section 14.11 This integrator is multisymplectic, 
i.e. the discrete multisymplectic form formula (I4.10P is satisfied. 

Example: Midpoint rule. In (|2.17p consider the 1-stage symplectic partitioned Runge-Kutta 
method with the coefficients an = a\\ = c\ = 1/2 and b\ = h\ = 1. This method is often called 
the midpoint rule and is a 2-nd order member of the Gauss family of quadratures. It can be easily 
shown (see |20| . |34| ) that the discrete Lagrangian (12.12P for this method is given by 

L d ( tj ,y>,t 1+1 , yi +l ) = At-L N ( » ±A , V _ I jt . + _ A tj , (4.18) 
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where At = tj+i — tj and y 3 = (yj, . . . , yii)- Using (|2,2p and (|2. lOj) we can write 



N 



L d (t ] ,yJ,t ]+1 ,y^) = Y,L(yi,yl +v yi:iy! +1 ), (4.19) 



8=0 



where we defined the discrete Lagrangian L: J l Y — > R by the formula 



L (yi^l + vyi+lyi +1 ) = Ai / c[(p(x),^ x (x),<p t (x),x,t j + ^At)dx dm 

with 



^ {x) = 2^x— + 2 Ax ' 

= W At Wi + i+ \ t l+1 Vi + i(*)- (4.21) 

Given the Lagrangian density C as in (|2.3[) . and assuming X(t,x) is known, one can evaluate 
the integral in (I4.20P explicitly It is now a straightforward calculation to show that the discrete 
variational principle ([4.17p for the discrete Lagrangian L as defined is equivalent to the Discrete 
Euler-Lagrange equations (|4.3p for L^, and consequently to (I2.17p . 

This shows that the 2-nd order Gauss method applied to (|2.17p defines a multisymplectic method 
in the sense of formula (|4,10p . However, for other symplectic partitioned Runge-Kutta methods of 
interest to us, namely the 4-th order Gauss and the 2-nd/4-th order Lobatto IIIA-IIIB methods, it 
is not possible to isolate a discrete Lagrangian L that would only depend on four values yj, 
y{i i j y{ + ■ The mentioned methods have more internal stages, and the equations ([2.17P couple 
them in a nontrivial way. Effectively, at any given time step the internal stages depend on all the 
values y{, . . . , yir and y{ +1 , . . . , y^ 1 , and it it not possible to express the discrete Lagrangian (|2.12p 
as a sum similar to (I4.19p . The resulting integrators are still variational, since they are derived by 
applying the discrete variational principle (|4.17p to some discrete action S, but this action cannot be 
expressed as the sum of L over all rectangles. Therefore, these integrators are not multisymplectic, 
at least not in the sense of formula (|4.10p . 



Constraints. Let the additional bundle be B = X x [0,X maa; ] and denote by X™ the elements of 
the fiber Bp. Define Y = Y (£> B. We have J k Y ^ J k Y (£> J k B. Suppose G : J k Y — ► E represents 
a discretization of the continuous constraint. For instance, one can enforce a uniform mesh by 
defining G : J l Y — > R, G{j l (p, j 1 X) = X x — 1 at the continuous level. The discrete counterpart will 
be defined on the discrete jet bundle J X Y by the formula 

G^,yl +1 ,yUl,yi + \xixi +1 ,xf+l,xi +1 ) = Xl+1 ~ Xl - 1. (4.22) 
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Arc-length equidistribution can be realized by enforcing ([2~24"p . that is, G : J$Y -> R, G(j$<p,j$X) = 
u 2 (p x ip xx + X^X^. The discrete counterpart will be defined on the discrete subbundle JqY by the 
formula 

G(y m i,X m r) = a 2 (y m t - y m2 ) 2 + (X m a - X m2 ) 2 - a 2 (y m 2 - y m i) 2 - (X m 2 - X m i) 2 , (4.23) 

where for convenience we used the notation introduced in (|4. 13j) and l,r = 1, . . . , 6. Note that A4.23P 
coincides with (j2.25|) . In fact, gi in ([2.25ft is nothing else but G computed on an element of JqY 
over the base 6-tuple □ such that Q] 2 = The only difference is that in (|2.25p we assumed gi 

might depend on all the field values at a given time step, while G only takes arguments locally, i.e. 
it depends on at most 6 field values on a given 6-tuple. 

A numerical scheme is now obtained by simultaneously solving the discrete Euler-Lagrange 
equations (|4.9p resulting from (|4.17p and the equation G = 0. If we know J/f -1 , Xj~ , yj and Xj 
for i = 1, . . . , N, this system of equations allows us to solve for ]£ + , X^ . This numerical scheme 
is multisymplectic in the sense similar to Proposition 12.41 If we take X(t,x) to be a sufficiently 
smooth interpolation of the values X? and substitute it in the problem (|4,17p . then the resulting 
multisymplectic integrator will yield the same numerical values yf + . 

4.3 Analysis of the Lagrange multiplier approach 
Continuous setting 

We now turn to describing the Lagrange multiplier approach in a multisymplectic setting. Similarly 
as in Section I4.2[ let the computational spacetime be X = R x [0, X max ] with coordinates (t, x) 
and consider the trivial configuration bundles tt^y '■ Y = X x R — > X and tt^b '■ B = X x 
[0,X max ] — > X. Let our scalar field be represented by a section <p : X — > Y with the coordinate 
representation <p(t, x) = (t, x, (p(t, x)) and our diffeomorphism by a section X : X — > B with the 
local representation X(t,x) = (t,x,X(t,x)). Let the total configuration bundle be Y = Y © B. 
Then the Lagrangian density (12.3|) can be viewed as a mapping C : J l Y = J^Y © J^B — > R. The 
corresponding action (jl.3p can now be expressed as 

S[<f,X) = [ C(j 1 (p,j 1 X)dtAdx, (4.24) 
Ju 

where U = [0, ^rnaa;] x [0, X max \. As before, the MMPDE constraint can be represented by a function 
G : J k Y — > R. Two sections (p and X satisfy the constraint if 

G(j k <p,j k X) = Q. (4.25) 

Vakonomic formulation. We now face the problem of finding the right equations of motion. 
We want to extremize the action functional (|4.24|) in some sense, subject to the constraint (|4.25p . 
Note that the constraint is essentially nonholonomic, as it depends on the derivatives of the fields. 
Assuming G is a submersion, G = defines a submanifold of J k Y, but this submanifold will not in 
general be the k-th jet of any subbundle of Y. Two distinct approaches are possible here. One could 
follow the Lagrange- d'Alembert principle and take variations of S first, but choosing variations V 
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(vertical vector fields on Y) such that the jet prolongations j k V are tangent to the submanifold 
G = 0, and then enforce the constraint G = 0. On the other hand, one could consider the variational 
nonholonomic problem (also called vakonomic), and minimize S over the set of all sections (ip,X) 
that satisfy the constraint G = 0, that is, enforce the constraint before taking the variations. If 
the constraint is holonomic, both approaches yield the same equations of motion. However, if the 
constraint is nonholonomic, the resulting equations are in general different. Which equations are 
correct is really a matter of experimental verification. It has been established that the Lagrange- 
d'Alembert principle gives the right equations of motion for nonholonomic mechanical systems, 
whereas the vakonomic setting is appropriate for optimal control problems (see [3], [4], [5]). 

We are going to argue that the vakonomic approach is the right one in our case. In Proposition [XT] 
we showed that in the unconstrained case extremizing S[<f)] with respect to <p was equivalent to 
extremizing S[<p, X] with respect to <p, and in Proposition 13.21 we showed that extremizing with 
respect to X did not yield new information. This is because there was no restriction on the fields (p 
and X, and for any given X there was a one-to-one correspondence between <p and <p given by the 
formula <p(t, x) = <p(t, X(t, x)), so extremizing over all possible (p was equivalent to extremizing over 
all possible <p. Now, let M be the set of all smooth sections (<p, X) that satisfy the constraint (14.25P 
such that X(t, .) is a diffeomorphism for all t. It should be intuitively clear that under appropriate 
assumptions on the mesh density function p, for any given smooth function <p(t,X), equation (I2.22P 
together with p(t,x) = (j)(t,X(t,x)) define a unique pair (tp,X) £ J\f (since our main purpose here 
is to only justify the application of the vakonomic approach, we do not attempt to specify those 
analytic assumptions precisely). Conversely, any given pair (<p,X) € N defines a unique function (j) 
through the formula cp(t,X) = (p(t,£(t,X)), where £(i, .) = X(t, .) , as in Section [3.11 Given this 
one-to-one correspondence and the fact that S[<f)] = S[<p,X] by definition, we see that extremizing S 
with respect to all smooth <fi is equivalent to extremizing S over all smooth sections (<p, X) € M. We 
conclude that the vakonomic approach is appropriate in our case, since it follows from Hamilton's 
principle for the original, physically meaningful, action functional S. 

Let us also note that our constraint depends on spatial derivatives only. Therefore, in the 
setting presented in Section [2] and Section [3] it can be considered holonomic, as it restricts the 
infinite-dimensional configuration manifold of fields that we used as our configuration space. In 
that case it is valid to use Hamilton's principle and minimize the action functional over the set of 
all allowable fields, i.e. those that satisfy the constraint G = 0. We did that by considering the 
augmented instantaneous Lagrangian (13.50p . 

In order to minimize S over the set of sections satisfying the constraint (|4.25p we will use the 
bundle-theoretic version of the Lagrange multiplier theorem, which we cite below after |32| . 

Theorem 4.1 (Lagrange multiplier theorem). Let ttm,£ '■ £ — ► -M. be an inner product bundle 
over a smooth manifold M., ^ a smooth section of tvm,£> an d h ■ M — > R a smooth function. 
Setting N = ^I /_1 (0) ; the following are equivalent: 

1. a G M is an extremum ofh\j\f, 

2. there exists an extremum a £ £ ofh:£ — > R such that ■km,e{&) = °~> 
where h(a) = h(ir Mt£ (a)) - {a^{-K M ,e{^))) e - 

Let us briefly review the ideas presented in |32j, adjusting the notation to our problem and 
generalizing when necessary. Let 
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C5° On = W = {<f, X) : U C X — > y } (4.26) 

be the set of smooth sections of TTyy on W. Then S : Cy(Y) — > R can be identified with /i in 
Theorem 14. 1| where .M = C^P(y). Furthermore, define the trivial bundle 

tt xv :V = X xR — > X (4.27) 

and let C^P(V) be the set of smooth sections A : U — > V, which represent our Lagrange multipliers 
and in local coordinates have the representation \(t,x) = (t,x, X(t,x)). The set C^P(V) is an inner 
product space with (Ai, A2) = fy A1A2 dt A dx. Take 

£ = C%(Y)xC%(V). (4.28) 
This is an inner product bundle over C^iY) with the inner product defined by 

((ff,Ai),((7,A 2 )) g = (Ai ) A 2 ). (4.29) 

We now have to construct a smooth section ^ : C^iY) — > E that will realize our constraint (14.25|) . 
Define the fiber- preserving mapping G : J k Y — > V such that for 1? G J k Y 

G(V)={ir XJk y(#),G(#))- (4-30) 

For instance, for k = 1, in local coordinates we have G(t,x,y,vt,v x ) = (t, x, G(t, x, y, vt , v x )). Then 
we can define 



*(a) = (a,Goj k a). (4.31) 

The set of allowable sections N C Cy(Y) is now defined by N = ^ _1 (0). That is, (<p,X) € N 
provided that G(j k ip,j k X) = 0. 

The augmented action functional 5c : £ — > R is now given by 



or denoting a = (jp, X, A) 



Sc[o\ = S[7r M , £ (a)} - (a, *(K M A*)))e> ( 4 ' 32 ) 



S c [<p, X, A] = S[<p, X] - (A, G o (j k Cp, j k X)) 



C{j l (p,j 1 X) dtAdx 



u 



X(t,x)G(j k ip,j k X)dt Adx 



u 



u 



£(j 1 ip,j 1 X)-\(t,x)G(j k <p,j k X) 



dt A dx. 



(4.33) 



Theorem 14.11 states, that if (<p,X, A) is an extremum of §c, then (jp,X) extremizes S over the set 
M of sections satisfying the constraint G = 0. Note that using the multisymplectic formalism we 
obtained the same result as (|3.5U|) in the instantaneous formulation, where we could treat G as a 
holonomic constraint. The dynamics is obtained by solving for a triple (if, X, A) such that 
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S c [vy o^^oI,^oA] = (4.34) 

=o 



for all r}y, rj^, r]y that keep the boundary conditions on dU fixed, where rf denotes the flow of 
vertical vector fields on respective bundles. 

Note that we can define Y c = Y © B © V and C c : J k % — > R by setting C c = C - X ■ G, 
i.e., we can consider a fc-th order field theory. If k = 1,2 then an appropriate multisymplectic form 
formula in terms of the fields ip, X and A will hold. Presumably, this can be generalized for k > 2 
using the techniques put forth in |28| . However, it is an interesting question whether there exists 
any multisymplectic form formula defined in terms of <p, X and objects on J k Y only. It appears 
to be an open problem. This would be the multisymplectic analog of the fact that the flow of 
a constrained mechanical system is symplectic on the constraint submanifold of the configuration 
space. 

Discretization 

Let us use the same discretization as discussed in Section I4~2l Assume we have a discrete Lagrangian 
L : J l Y — > R, the corresponding discrete action S[<^,X], and a discrete constraint G : J l Y — > R 
or G : JqY — > R. Note that S is essentially a function of 2MN variables and we want to extremize 
it subject to the set of algebraic constraints G = 0. The standard Lagrange multiplier theorem 
proved in basic calculus textbooks applies here. However, let us work out a discrete counterpart of 
the formalism introduced at the continuous level. This will facilitate the discussion of the discrete 
notion of multisymplecticity. Let 

C U (Y) = {a = (<p, X) : U C X — > Y} (4.35) 

be the set of discrete sections of 7T^y : Y — > X. Similarly, define the discrete bundle V = X xR and 
let Cu (V) be the set of discrete sections A : IAq — > V representing the Lagrange multipliers, where 
C U is defined below. Let = (J,i, X(J,i)) with A] = be the local representation. 

The set Cu (V) is an inner product space with (A, p) = J2(ji)<=u HfA- Take £ = Cu(Y) x Cu Q (V). 
Just like at the continuous level, £ is an inner product bundle. However, at the discrete level it is 
more convenient to define the inner product on £ in a slightly modified way. Since there are some 
nuances in the notation, let us consider the cases k = 1 and k = 2 separately. 

Case k=l. Let U = {(j, i) eU\j < M, i < N}. Define the trivial bundle V = X u x R and let 
C u n{V) be the set of all sections of V defined on IA U . For a given section A S Cu {V) we define its 
extension A G C u n (V) by 

A(D) = (n i A(D 1 )), (4.36) 

that is, A assigns to the square □ the value that A takes on the first vertex of that square. Note 
that this operation is invertible: given a section of C u n (V) we can uniquely determine a section of 
Cua(y)- AVe can define the inner product 

(A,A) = ^A^KD 1 ). (4.37) 
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One can easily see that we have (A, ft) = (A, p), so by a slight abuse of notation we can use the same 
symbol (.,.) for both inner products. It will be clear from the context which definition should be 
invoked. We can now define an inner product on the fibers of £ as 

((a,\),(a,jl)) £ = (X,fi} = (\,ii}. (4.38) 

Let us now construct a section ^ : Cu(Y) — > £ that will realize our discrete constraint G. First, 
in analogy to (14.301) . define the fiber- preserving mapping G : J^Y — > V such that 

G(y nl ,X D r) = (D,G(y Dh X a r)), (4.39) 

where l,r = 1,2,3,4. We now define ^ by requiring that for a 6 Cu(Y) the extension (|4.36p of 
^f(cr) is given by 

#(ct) = {a,Goj l a). (4.40) 

The set of allowable sections J\f C CuiY) is now defined by N = ^ _1 (0)— that is, (<p,X) £ N 
provided that G^ 1 ^, ^X) = for all □ € U u . The augmented discrete action Sq '■ £ — > K is 
therefore 

Sc[v,\] = S[<r]-((<T,\)M<r)) e 
= S[<j]-(\,Goj 1 a) 

= y, lu 1 *) - E A ( nl MA) 

UCU UcU 

= E (lO'VJ-A^JGO'V)). (4.41) 

UcU 



By the standard Lagrange multiplier theorem, if (<f,X, A) is an extremum of Sc, then (cp,X) is an 
extremum of S over the set Af of sections satisfying the constraint G = 0. The discrete Hamilton 
principle can be expressed as 

S C [<Pe,XeX]=0 (4-42) 

for all vector fields V on Y, W on B, and Z on V that keep the boundary conditions on dU fixed, 
where <p e (J,i) = F^ i% (<p(j,i)) and F^ 1 is the flow of Vji on M, and similarly for X e and A e . The 
discrete Euler-Lagrange equations can be conveniently computed if in (|4.42|) one focuses on some 
£ intZY. With the convention <p(j,i) = yj, X(j,i) = X\. A(j, i) = \\, we write the terms of 
Sc containing y\ . X\ and \\ explicitly as 



d_ 

Ik 
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S c = ...+L{ylyi +v y{^,yt\xi,Xi +1 ,Xi^,Xi +1 ) 

+ L{^ 1 M,yt\y{-lxi 1 ,xi,xi +1 ,xi^) 

• 1(4 ! . !f j. !f j ,..v/ /..v/ ; ..v/..v/ ,) 

• M/// ,.///• .V;' l ..V/. ; ; ..V/. ; ..V/) 

• A^( / ^ // f. l ., / /;;. // r l ..v/..v/.,..v/. l l ..v/ 1 ) 

• A' v MM + \v£lM ; ..V/..V/- | ..V;' ,') 

. V;( y/ / L \yf '.^ ,..v/ ; ; ..v/ -..v/..v/ ,) 

• a^ ; ..v/. ; ; ..v/.,..v/) • ... 



(4.43) 



The discrete Euler-Lagrange equations are obtained by differentiating with respect to y\, X\ and 
A^ , and can be written compactly as 



E 

l,D 

UA=a l 



dL 



dG 
dy l 



E 

i,a 



dL 
dX l 



dG 



+ ^□ 1 ^r(yn 1 ' • • • 'yn 4 )^ 1 ! • • • )^n 4 ) 



0, 



G(jjlyl +1 ,yi:iyl +1 ,xi,xi +1 ,x£,xr)=o 



1.44) 



for all (j, i) G intW. If we know y| 1 , _X| , yj, Xf and A^ 1 for i = 1, ...,N, this system of 
equations allows us to solve for yf + , Xf +1 and Xj. 

Note that we can define Yo = Y © ,8 © V and the augmented Lagrangian Lc : J 1 ^* — ^ ^ by 
setting 



(4.45) 



that is, we can consider an unconstrained field theory in terms of the fields <p, X and A. Then, the 
solutions of (|4.44p satisfy the multisymplectic form formula (14.101) in terms of objects defined on 
J l Y C - 



Case k=2. Let Uq = {(j, i) G U \ j < M, 1 < i < N}. Define the trivial bundle V = X m x R and 
let CymiV) be the set of all sections of V defined on £Y m . For a given section A G C^ (V) we define 
its extension A G C w m (V) by 
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A(m) = (m,A(m 2 )), 



(4.46) 



that is, A assigns to the 6-tuple □ the value that A takes on the second vertex of that 6-tuple. Like 
before, this operation is invertible. We can define the inner product 

(A,A)= £ A(m 2 V(m 2 ) (4.47) 

mcu 

and the inner product on £ as in (|4.38p . Define the fiber- preserving mapping G : JqY — > V such 
that 

G(y ml ,X mr ) = (m,G(y ml ,X m r)), (4.48) 

where I, r = 1, . . . , 6. We now define ^ by requiring that for a 6 Cu(Y) the extension (|4.46p of ^(c) 
is given by 

4>(a) = (a,Gof a). (4.49) 

Again, the set of allowable sections is N = ^ ,_1 (0). That is, (</>, X) € M provided that G(j 2 ^, 3qX) = 
for all □ E W m . The augmented discrete action So : £ — > M is therefore 

Sbk,A] = S>]-((a,A),*(a))^ 
= 5H-(A,Goj 2 (T ) 

= £ L(iV) - £ A(m 2 )G(j M. (4.50) 

□cw mew 

Writing out the terms involving yj, Xj and A| explicitly, as in (|4,43p . and invoking the discrete 
Hamilton principle (|4.42p . one obtains the discrete Euler-Lagrange equations, which can be com- 
pactly expressed as 
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dL 

^ -jr-iivu 1 -, ■ ■ ■ lUnfiiXtii, . . . ,X n *)+ 



dy l 



dG 



EdL 



W,i)=n ! 



</(,,/ ,-///-///. : . . yy/ '-/// - f - - V / , . .V/. .V/. ; . .V// ; : . .V/ ' ; . .V/ ' ) (4.51) 

for all G intZY. If we know , -X|' — , yf, Xj and A] 1 for « = 1,...,N, this system of 

equations allows us to solve for yj + , and Aj. 

Let us define the extension L ex t : </q^ — ^ ^ °^ the Lagrangian density L by setting 



{ L(y nl ,...,X D ,) if m 2 = (j,0),(i,iV + l), 

where n = mn^, (4.52) 
, lEncm^n 1 '--- > X n 4 ) otherwise. 



Let us also set G(y D i,. . . ,X n4 ) = if m 2 = (j,0), (j, iV + 1). Define A = {m | m 2 ,m 5 G ^/}. Then 
(|4,50p can be written as 



S c [a, A] = E [ Z ext(io 2 ^) " A(m 2 )G(j M] = LctifajlX), (4.53) 
meA meA 

where the last equality defines the augmented Lagrangian Lc '■ JqYc — > K for Yq = Y © B © V. 
Therefore, we can consider an unconstrained second-order field theory in terms of the fields (p, X 
and A, and the solutions of (|4,5ip will satisfy a discrete multisymplectic form formula very similar 
to the one proved in [28J. The only difference is the fact that the authors analyzed a discretization 
of the Camassa-Holm equation and were able to consider an even smaller subbundle of the second 
jet of the configuration bundle. As a result it was sufficient for them to consider a discretization 
based on squares □ rather than 6-tuples □. In our case there will be six discrete 2- forms for 
I = 1, ... ,6 instead of just four. 

Remark. In both cases we showed that our discretization leads to integrators that are multisym- 
plectic on the augmented jets J k Yc- However, just like in the continuous setting, it is an interesting 
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problem whether there exists a discrete multisymplectic form formula in terms of objects defined 
on J k Y only. 



Example: Trapezoidal rule. Consider the semi-discrete Lagrangian (|3.8p . We can use the 

trapezoidal rule to define the discrete Lagrangian (|2.1ip as 

( , 54) 

where y J = (t/J, . . . , y J N ) and X d = (X{, . . . , X° N ). The constrained version (see |34| ) of the Discrete 
Euler-Lagrange equations (|4.3p takes the form 



D 2 L d (qi-\q3) + D x L d {q> ,qi +1 ) = Dg(<?) T \i, 

9(q 1+1 ) = 0, (4.55) 

where for brevity q 3 = (yj, X^, . . . , y 3 N , X 3 N ), X J = (A{, . . . , Ajy) and g is an adaptation constraint, 
for instance (|2.25p . If qi~ , qi are known, then (|4,55p can be used to compute q 3+ and A J . It is 
easy to verify that the condition (I3.58P is enough to ensure solvability of (14. 55|) . assuming the time 
step At is sufficiently small, so there is no need to introduce slack degrees of freedom as in (|3.59p . 
If the mass matrix (I3.1ip was constant and nonsingular, then (|4,55p would result in the SHAKE 
algorithm, or in the RATTLE algorithm if one passes to the position-momentum formulation (see 
[20], M)- 

Using (13. 2p and (13. 5p we can write 

N 

L d ( y'. X'. y*+\ A- • ; ; .V/. .V/. , . .V/.' . .V/ ' ; ) . (4.56) 

i=0 

where we defined the discrete Lagrangian L : J X Y — > R by the formula 



f (iP iP iP +l u j+1 X j X j X j+1 X j+1 
At f Xi + l 



2 



£ (f 3 (x) , X J (x) , (pi (x) , XI (x) , ip t (x) , X t (x) )dx 



+ T J £(^ + Hx),Xi + \x),rt + H X ),XtH X ),<p t (x),X t (x))dx (4.57) 
with 
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ifP{x) = y{r}i(x) +y{ +1 r)i +1 (x), 
j j 

-j , \ Vi+l ~ Vi 

aw = — s^— • 

J+ 1 _ J J. +1 - J 

<Pt(x) = Vi At + i+ \ t l+1 Vi + i(x), (4.58) 

and similarly for X(x). Given the Lagrangian density C as in fl3.Tf> one can compute the integrals 
in (|4.57p explicitly. Suppose that the adaptation constraint g has a 'local' structure, for instance 

gttf,*) = r;^. // /. ; ., / f; l l . // f ; ..v/..v/. l ..v/. l l ..v/- ; ). (4.59) 

as in (|4^2j) or 

gi (y j ,Xi) = G(y mh X m r), where m 2 = (4.60) 

as in (I4.23p . It is straightforward to show that (|4.44j) or (|4.5ip are equivalent to (|4.55p . that is, the 
variational integrator defined by (I4.55P is also multisymplectic. 

For reasons similar to the ones pointed out in Section 14. 2| the 2-nd and 4-th order Lobatto 
IIIA-IIIB methods that we used for our numerical computations are not multisymplectic. 



5 Numerical results 

5.1 The Sine-Gordon equation 

We applied the methods discussed in the previous sections to the Sine-Gordon equation 

d 2 4> d 2 4> , , 

^-^+ w =°- (5 ' 1} 

This equation results from the (l+l)-dimensional scalar field theory with the Lagrangian density 

C(ct>AxAt) = \$ - \<t>x " (1 " cos^). (5.2) 

The Sine-Gordon equation arises in many physical applications. For instance, it governs the prop- 
agation of dislocations in crystals, the evolution of magnetic flux in a long Josephson-j unction 
transmission line or the modulation of a weakly unstable baroclinic wave packet in a two-layer fluid. 
It also has applications in the description of one-dimensional organic conductors, one-dimensional 
ferromagnets, liquid crystals, or in particle physics as a model for baryons (see |45|). 

The Sine-Gordon equation has interesting soliton solutions. A single soliton traveling at the 
speed v is given by 



4>s(X,t) = 4arctan 



exp 



X-Xo-vt 



(5.3) 
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Figure 5.1: The single soliton solution of the Sine-Gordon equation. 

It is depicted in Figure [5TT1 The backscattering of two solitons, each traveling with the velocity v, 
is described by the formula 



cj>Ss(X,t) = 4arctan 



v sinh( 



x 



cosh( 



vt 



VT- 



(5.4) 



It is depicted in Figure 15.21 Note that if we restrict X > 0, then this formula also gives a single 
soliton solution satisfying the boundary condition 0(0, t) = 0, that is, a soliton bouncing from a 
rigid wall. 

5.2 Generating consistent initial conditions 

Suppose we specify the following initial conditions 



<KX,0) = a(X), 
<Pt(X,0)=b(X), 



(5.5) 



and assume they are consistent with the boundary conditions (jl.2p . In order to determine appro- 
priate consistent initial conditions for (|2.15p and (|3.62|) we need to solve several equations. First 



we solve for the y^s and Xj's. We have yo 
determined by solving the system 



^Li Vn+i = <t>Ri Xq = 0, Xn + \ = X„ 



The rest are 



yi = a(Xi), 

= • • • ,Vn,Xi, . . . ,X N ), 



(5.6) 
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Figure 5.2: The two-soliton solution of the Sine-Gordon equation. 
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for i = 1,...,N. This is a system of 2N nonlinear equations for 2N unknowns. We solve it 
using Newton's method. Note, however, that we do not a priori know good starting points for 
Newton's iterations. If our initial guesses are not close enough to the desired solution, the iterations 
may converge to the wrong solution or may not converge at all. In our computations we used the 
constraints (I2.25P . We found that a very simple variant of a homotopy continuation method worked 
very well in our case. Note that for a = the set of constraints (|2,25p generates a uniform mesh. 
In order to solve (15. 6ft for some a > 0, we split [0, a] into d subintervals by picking = (k/d) - a 
for k = l,...,d. We then solved fj5.6|) with oi\ using the uniformly spaced mesh points xj® = 
(i/(N + 1)) • X max as our initial guess, resulting in X- and . Then we solved (|5.6|) with o?2 
using X^ and y^' as the initial guesses, resulting in xf^ and y\ . Continuing in this fashion, 
we got xj^ and y^ as the numerical solution to (|5.6|) for the original value of a. Note that for 
more complicated initial conditions and constraint functions, predictor-corrector methods should be 
used — see pQ for more information. Another approach to solving (I5.6P could be based on relaxation 
methods (see [7|, [25]). 

Next, we solve for the initial values of the velocities y% and X{. Since ip(x,t) = 4>(X(x, t), t), 
we have <pt(x,t) = 4>x{X(x,t),t)Xt(x,t) + (f>t(X(x, t), t). We also require that the velocities be 
consistent with the constraints. Hence the linear system 



yi = a'(Xi)Xi + b(Xi), i = l,...,N 

= ^(y,X)y+gL(y,X)X. (5.7) 

This is a system of 2N linear equations for the 2N unknowns y^ and Xi, where y = (yi,. . . , J/jv) 
and X = (Xi, . . . ,Xjy). We can use those velocities to compute the initial values of the conjugate 
momenta. For the control-theoretic approach we use pi = dL^/dyi, as in Section 12. 3( and for 
the Lagrange multiplier approach we use (|3,10p . In addition, for the Lagrange multiplier approach 
we also have the initial values for the slack variables = and their conjugate momenta Bi = 
dL^jdfi = 0. It is also useful to use (I3.57P to compute the initial values of the Lagrange multipliers 
Aj that can be used as initial guesses in the first iteration of the Lobatto IIIA-IIIB algorithm. The 
initial guesses for the slack Lagrange multipliers are trivially /ij = 0. 



5.3 Convergence 

In order to test the convergence of our methods as the number of mesh points is increased, we 
considered a single soliton bouncing from two rigid walls at X = and X = X max = 25. We 
imposed the boundary conditions 0i = and 0# = 2ir, and as initial conditions we used (|5.3p with 
Xq = 12.5 and v = 0.9. It is possible to obtain the exact solution to this problem by considering 
a multi-soliton solution to (|5.ip on the whole real line. Such a solution can be obtained using a 
Backlund transformation (see |llj . |45j). However, the formulas quickly become complicated and, 
technically, one would have to consider an infinite number of solitons. Instead, we constructed a 
nearly exact solution by approximating the boundary interactions with <\5A\i : 



'exact 



(X,t) 



<j>ss{X- 
<j>ss(X,t 



X m ax j t 



(An + 1)T) + 2vr if t € 
(4n + 3)T) if t G 



4nT, (4n + 2)T) , 

(4n + 2)T, (4n + 4)T), 



(5.8) 
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where n is an integer number and T satisfies 4>ss(X ma x/2, T) = tt (we numerically found T ~ 13.84). 
Given how fast (|5,3p and f|5.4[) approach its asymptotic values, one may check that (|5.8p can be 
considered exact to machine precision. 

We performed numerical integration with the constant time step At = 0.01 up to the time 
T m ax = 50. For the control-theoretic strategy we used the 1-stage and 2-stage Gauss method 
(2-nd and 4-th order respectively), and the 2-stage and 3-stage Lobatto IIIA-IIIB method (also 
2-nd/4-th order). For the Lagrange multiplier strategy we used the 2-stage and 3-stage Lobatto 
IIIA-IIIB method for constrained mechanical systems (2-nd/4-th order). See |20| . |21| . |22| for 
more information about the mentioned symplectic Runge-Kutta methods. We used the constraints 
(|2.25|) based on the generalized arclength density (I2.23p . We chose the scaling parameter to be 
a = 2.5, so that approximately half of the available mesh points were concentrated in the area of 
high gradient. A few example solutions are presented in Figure I5.3H5.6I Note that the Lagrange 
multiplier strategy was able to accurately capture the motion of the soliton with merely 17 mesh 
points (that is, N = 15). The trajectories of the mesh points for several simulations are depicted 
in Figure 15.81 and Figure 15.91 An example solution computed on a uniform mesh is depicted in 
Figure [5771 

For the convergence test, we performed simulations for several N in the range 15-127. For 
comparison, we also computed solutions on a uniform mesh for N in the range 15-361. The numerical 
solutions were compared against the solution (15. 8j) . The L°° errors are depicted in Figure [5.101 The 
L°° norms were evaluated over all nodes and over all time steps. Note that in case of a uniform mesh 
the spacing between the nodes is Ax = X max / (A+l), therefore the errors are plotted versus (A+l). 
The Lagrange multiplier strategy proved to be more accurate than the control-theoretic strategy. As 
the number of mesh points is increased, the uniform mesh solution becomes quadratically convergent, 
as expected, since we used linear finite elements for spatial discretization. The control-theoretic 
strategy also shows near quadratic convergence, whereas the Lagrange multiplier method seems to 
converge slightly slower. While there are very few analytical results regarding the convergence of 
r-adaptive methods, it has been observed that the rate of convergence depends on several factors, 
including the chosen mesh density function. Our results are consistent with the convergence rates 
reported in [2j and |50| . Both papers deal with the viscous Burgers' equation, but consider different 
initial conditions. Computations with the arclength density function converged only linearly in |2j, 
but quadratically in [50j. 




5.4 Energy conservation 

As we pointed out in Section 12. 6[ the true power of variational and symplectic integrators for 
mechanical systems lies in their excellent conservation of energy and other integrals of motion, even 
when a big time step is used. In order to test the energy behavior of our methods, we performed 
simulations of the Sine-Gordon equation over longer time intervals. We considered two solitons 
bouncing from each other and from two rigid walls at X = and X max = 25. We imposed 
the boundary conditions <pL = —2ir and <j)R = 2n, and as initial conditions we used 4>(X, 0) = 
(pSs(X — 12.5, —5) with v = 0.9. We ran our computations on a mesh consisting of 27 nodes 
(N=25). Integration was performed with the time step At = 0.05, which is rather large for this 
type of simulations. The scaling parameter in (|2.25p was set to a = 1.5, so that approximately half 
of the available mesh points were concentrated in the areas of high gradient. An example solution 
is presented in Figure 15.111 

The exact energy of the two-soliton solution can be computed using ([2.4p . It is possible to 
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Figure 5.3: The single soliton solution obtained with the Lagrange multiplier strategy for N = 15. 
Integration in time was performed using the 4-th order Lobatto IIIA-IIIB scheme for constrained 
mechanical systems. The soliton moves to the right with the initial velocity v = 0.9, bounces from 
the right wall at t = 13.84 and starts moving to the left with the velocity v = —0.9, towards the 
left wall, from which it bounces at t = 41.52. 
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Figure 5.4: The single soliton solution obtained with the Lagrange multiplier strategy for N = 31. 
Integration in time was performed using the 4-th order Lobatto IIIA-IIIB scheme for constrained 
mechanical systems. 



50 





Figure 5.5: The single soliton solution obtained with the control-theoretic strategy for N = 22. 
Integration in time was performed using the 4-th order Gauss scheme. Integration with the 4-th 
order Lobatto IIIA-IIIB yields a very similar level of accuracy. 
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Figure 5.6: The single soliton solution obtained with the control-theoretic strategy for N = 31. 
Integration in time was performed using the 4-th order Gauss scheme. Integration with the 4-th 
order Lobatto IIIA-IIIB yields a very similar level of accuracy. 
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Figure 5.7: The single soliton solution computed on a uniform mesh with N = 31. Integration in 
time was performed using the 4-th order Gauss scheme. Integration with the 4-th order Lobatto 
IIIA-IIIB yields a very similar level of accuracy. 
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Figure 5.8: The mesh point trajectories (with zoomed-in insets) for the Lagrange multiplier strategy 
for N = 22 (left) and N = 31 (right). Integration in time was performed using the 4-th order Lobatto 
IIIA-IIIB scheme for constrained mechanical systems. 




Figure 5.9: The mesh point trajectories (with zoomed-in insets) for the control-theoretic strategy 
for N = 22 (left) and N = 31 (right). Integration in time was performed using the 4-th order Gauss 
scheme. Integration with the 4-th order Lobatto IIIA-IIIB yields a very similar result. 
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Figure 5.10: Comparison of the convergence rates of the discussed methods. Integration in time 
was performed using the 4-th order Lobatto IIIA-IIIB method for constrained systems in case of the 
Lagrange multiplier strategy, and the 4-th order Gauss scheme in case of both the control-theoretic 
strategy and the uniform mesh simulation. The 4-th order Lobatto IIIA-IIIB scheme for the control- 
theoretic strategy and the uniform mesh simulation yields a very similar level of accuracy. Also, 
using 2-nd order integrators gives very similar error plots. 
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Figure 5.11: The two-soliton solution obtained with the control-theoretic and Lagrange multiplier 
strategies for N = 25. Integration in time was performed using the 4-th order Gauss quadrature 
for the control-theoretic approach, and the 4-th order Lobatto IIIA-IIIB quadrature for constrained 
mechanical systems in case of the Lagrange multiplier approach. The solitons initially move towards 
each other with the velocities v = 0.9, then bounce off of each other at t = 5 and start moving 
towards the walls, from which they bounce at t = 18.79. The solitons bounce off of each other 
again at t = 32.57. This solution is periodic in time with the period T per j ^ = 27.57. The nearly 
exact solution was constructed in a similar fashion as (15. 8j) . As the simulation progresses, the 
Lagrange multiplier solution gets ahead of the exact solution, whereas the control-theoretic solution 
lags behind. 
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compute that integral explicitly to obtain E = 16/\/l — v 2 ~ 36.71. The energy associated with 
the semi-discrete Lagrangian (|3.8p can be expressed by the formula 



E N = -q T M N (q)q + R N {q), 



(5.9) 



where Rn was defined in (|3.52p and for our Sine-Gordon system is given by 



N r 



1 ( Vk+i - Vk 

2 V^fc+i ~~ Xj- 




rn(q) = 



+ i 



sin y k+ i - sin y k 



(^fc+i — x k ), 



(5.10) 



fc=0 L 



y k+1 - yk 



and Mjv is the mass matrix (I3.1ip . The energy En is an approximation to (|2.4I) if the integrand is 
sampled at the nodes Xq,. . .,Xn+i and then piecewise linearly approximated. Therefore, we used 
En to compute the energy of our numerical solutions. 

The energy plots for the Lagrange multiplier strategy are depicted in Figure 15.121 We can see 
that the energy stays nearly constant in the presented time interval, showing only mild oscillations, 
which are reduced as higher order of integration in time is used. The energy plots for the control- 
theoretic strategy are depicted in Figure 15.131 In this case the discrete energy is more erratic and 
not as nearly preserved. Moreover, the symplectic Gauss and Lobatto methods show virtually the 
same energy behavior as the non-symplectic Radau IIA method, which is known for its excellent 
stability properties when applied to stiff differential equations (see |22|). It seems that we do not 
gain much by performing symplectic integration in this case. It is consistent with our observations 
in Section 12.61 and shows that the control-theoretic strategy does not take the full advantage of the 
underlying geometry. 

As we did not use adaptive time-stepping and did not implement any mesh smoothing techniques, 
the quality of the mesh deteriorated with time in all the simulations, eventually leading to mesh 
crossing, i.e. two mesh points collapsing or crossing each other. The control-theoretic strategy, even 
though less accurate, retained good mesh quality longer, with the break-down time Tbreak 

> 1000, 

as opposed to Ti, reak ~ 600 in case of the Lagrange multiplier approach (both using a rather large 
constant time step). We discuss extensions to our approach for increased robustness in Section [6j 

6 Summary and future work 

We have proposed two general ideas how r-adaptive meshes can be applied in geometric numerical 
integration of Lagrangian partial differential equations. We have constructed several variational and 
multisymplectic integrators and discussed their properties. We have used the Sine-Gordon model 
and its solitonic solutions to test our integrators numerically. 

Our work can be extended in many directions. Interestingly, it also opens many questions 
in geometric mechanics and multisymplectic field theory. Addressing those questions will have a 
broader impact on the field of geometric numerical integration. 

Non-hyperbolic equations 

The special form of the Lagrangian density (13.6|) we considered leads to a hyperbolic PDE, which 
poses a challenge to r-adaptive methods, as at each time step the mesh is adapted globally in 
response to local changes in the solution. Causality and the structure of the characteristic lines 
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Figure 5.12: The discrete energy En for the Lagrange multiplier strategy. Integration in time was 
performed with the 2-nd (top) and 4-th (bottom) order Lobatto IIIA-IIIB method for constrained 
mechanical systems. The spikes correspond to the times when the solitons bounce off of each other 
or of the walls. 
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Figure 5.13: The discrete energy En for the control-theoretic strategy. Integration in time was per- 
formed with the 4-th order Gauss (top), 4-th order Lobatto IIIA-IIIB (middle) and non-symplectic 
5-th order Radau IIA (bottom) methods. 
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of hyperbolic systems make r-adaptation prone to instabilities and integration in time has to be 
performed carefully. The literature on r-adaptation almost entirely focuses on parabolic problems 
(see [7], [25j and references therein). Therefore, it would be interesting to apply our methods to 
PDEs that are first-order in time, for instance the Korteweg-de Vries, Nonlinear Schrodinger or 
Camassa-Holm equations. All three equations are first-order in time and are not hyperbolic in 
nature. Moreover, all can be derived as Lagrangian field theories (see [Sj, [9], |1Q| . [14] , |17| . 
|28j). The Nonlinear Schrodinger equation has applications to optics and water waves, whereas the 
Korteweg-de Vries and Camassa-Holm equations were introduced as models for waves in shallow 
water. All equations possess interesting solitonic solutions. The purpose of r-adaptation would be 
to improve resolution, for instance, to track the motion of solitons by placing more mesh points 
near their centers and making the mesh less dense in the asymptotically flat areas. 

Hamiltonian Field Theories 

Variational multisymplectic integrators for field theories have been developed in the Lagrangian 
setting (|28j, |31j). However, many interesting field theories are formulated in the Hamiltonian 
setting. They may not even possess a Lagrangian formulation. It would be interesting to con- 
struct Hamiltonian variational integrators for multisymplectic PDEs by generalizing the variational 
characterization of discrete Hamiltonian mechanics. This would allow to handle Hamiltonian PDEs 
without the need for converting them to the Lagrangian framework. Recently Leok & Zhang |29) 
and Vankerschaver & Ciao &: Leok |49j have laid foundations for such integrators. It would also 
be interesting to see if the techniques we used in our work could be applied in order to construct 
r-adaptive Hamiltonian integrators. 

Time adaptation based on local error estimates 

One of the challenges of r-adaptation is that it requires solving differential-algebraic or stiff ordinary 
differential equations. This is because there are two different time scales present: one defined by the 
physics of the problem and one following from the strategy we use to adapt the mesh. Stiff ODEs 
and DAEs are known to require time integration with an adaptive step size control based on local 
error estimates (see [6], |22j). In our work we used constant time-stepping, as adaptive step size 
control is difficult to combine with geometric numerical integration. Classical step size control is 
based on past information only, time symmetry is destroyed and with it the qualitative properties of 
the method. Hairer & Soderlind [24J developed explicit, reversible, symmetry-preserving, adaptive 
step size selection algorithms for geometric integrators, but their method is not based on local error 
estimation, thus it is not useful for r-adaptation. Symmetric error estimators are considered in |27j 
and some promising results are discussed. Hopefully, the ideas presented in those papers could be 
combined and generalized. The idea of Asynchronous Variational Integrators (see |30| ) could also 
be useful here, as this would allow to use a different time step for each cell of the mesh. 

Constrained multisymplectic field theories 

The multisymplectic form formula (|4.6p was first introduced in |31| . The authors, however, consider 
only unconstrained field theories. In our work we start with the unconstrained field theory (jl.lj) . 
but upon choosing an adaptation strategy represented by the constraint G = we obtain a con- 
strained theory, as described in Section [3] and Section 14.31 Moreover, this constraint is essentially 
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nonholonomic, as it contains derivatives of the fields, and the equations of motion are obtained 
using the vakonomic approach (also called variational nonholonomic) rather than the Lagrange- 
d'Alembert principle. All that gives rise to many very interesting and general questions. Is there a 
multisymplectic form formula for such theories? Is it derived in a similar fashion? Do variational 
integrators obtained this way satisfy some discrete multisymplectic form formula? These issues have 
been touched upon in [32], but by no means resolved. 

Mesh smoothing and variational nonholonomic integrators 

The major challenge of r-adaptive methods is mesh crossing, which occurs when two mesh points 
collapse or cross each other. In order to avoid mesh crossing and retain good mesh quality, mesh 
smoothing techniques were developed ([7], |25j). They essentially attempt to regularize the exact 
equidistribution constraint G = by replacing it with the condition edX/dt = G, where e is 
a small parameter. This can be interpreted as adding some attraction and repulsion pseudoforces 
between mesh points. If one applies the Lagrange multiplier approach to r-adaptation as described in 
Section [3J then upon finite element discretization one obtains a finite dimensional Lagrangian system 
with a nonholonomic constraint. This constraint is enforced using the vakonomic (nonholonomic 
variational) formulation. Variational integrators for systems with nonholonomic constraints have 
been developed mostly in the Lagrange-d'Alembert setting. Surprisingly, there seems to be virtually 
no literature regarding discrete vakonomic mechanics. In a recent paper, Garcia &: Fernandez Sz 
Rodrigo [15J address variational integrators for vakonomic systems. The ideas presented in that 
paper could be used to design structure-preserving mesh smoothing techniques. 
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